PET-enabled Dual-Energy CT: Image Reconstruction and A Proof-of-Concept Computer Simulation Study

As illustrated in Fig. 1, a standard PET/CT scan normally consists of a PET emission scan at 511 keV and a X-ray CT transmission scan commonly acquired at 80-140 kVp. The X-ray CT image has been mainly used for PET attenuation correction with which the PET scan provides a functional image describing the radiotracer distribution in the subject. Our proposed method exploits the potential of a standard PET emission scan for high-energy GCT imaging. Different from X-ray CT which uses an external X-ray source to generate tomographic data, here the PET-enabled GCT exploits the internal “$\gamma$-rays” generated by annihilation radiation of PET radiotracer decay in the object. The GCT image obtained from PET is then combined with the low-energy X-ray CT to form dual-energy CT imaging.

There are potentially multiple methods for obtaining a GCT image. In this work, we exploit a PET attenuation-activity joint reconstruction method, as described below.

PET projection measurement $\bm{y}$ can be well modeled as independent Poisson random variables with the log-likelihood function,

$$L(\bm{y}|\bm{\lambda},\bm{\mu})=\sum_{i=1}^{N_{d}}\sum_{m=1}^{N_{t}}y_{i,m}\log\bar{y}_{i,m}(\bm{\lambda},\bm{\mu})-\bar{y}_{i,m}(\bm{\lambda},\bm{\mu}),$$

(1)

where $i$ denotes the index of PET detector pair and $m$ denotes the index of time-of-flight (TOF) bin. $N_{d}$ is the total number of detector pairs and $N_{t}$ is the number of TOF bins. The expectation of the PET projection data $\bar{\bm{y}}(\bm{\lambda},\bm{\mu})$ is related to the activity image $\bm{\lambda}$ and attenuation image $\bm{\mu}$ at 511 keV through

$$\bar{\bm{y}}_{m}(\bm{\lambda},\bm{\mu})=\mathrm{diag}\left\{\bm{n}_{m}(\bm{\mu})\right\}\bm{G}_{m}\bm{\lambda}+\bm{r}_{m}$$

(2)

where $\bm{G}_{m}$ is the PET detection probability matrix for the $m$th timing bin and $\bm{r}$ accounts for the expectation of random and scattered events. $\bm{n}_{m}(\bm{\mu})$ is the normalization factor for TOF bin $m$, of which the $i$th element is

$$n_{i,m}(\bm{\mu})=c_{i,m}\cdot\exp{\big{(}-[\bm{A}\bm{\mu})]_{i}\big{)}}$$

(3)

where $c_{i,m}$ denotes the multiplicative factor other than the attenuation correction factor and $\bm{A}$ is the system matrix for transmission imaging.

For standard PET/CT imaging, the attenuation image $\bm{\mu}$ is normally predetermined from a X-ray CT scan and the PET reconstruction problem only estimates the $\bm{\lambda}$ image [37]. $\bm{\mu}$ can be approximated from a X-ray CT scan using a bilinear scaling conversion of linear attenuation coefficient from the X-ray energy (e.g., 140 kVp) to 511 keV [38].

The maximum-likelihood attenuation and activity (MLAA) estimation method [35, 9] seeks the estimates of both $\bm{\mu}$ and $\bm{\lambda}$ simultaneously by maximizing the Poisson log-likelihood,

$$\hat{\bm{\lambda}},\hat{\bm{\mu}}=\arg\max_{\bm{\lambda}\geq\bm{0},\bm{\mu}\geq\bm{0}}L(\bm{y}|\bm{\lambda},\bm{\mu}).$$

(4)

The MLAA formulation was first proposed for non-TOF data [35] but the simultaneous estimation suffers from cross-talk artifacts despite some encouraging results [35, 39]. The method was later demonstrated more effective for TOF data [40, 41]. A seminal theoretical work later proved that TOF data determine $\bm{\mu}$ up to a constant [8, 9]. Since then, the MLAA method has received a wide range of interests (e.g., [17, 12, 11, 21, 14, 19, 18, 13, 22, 23, 15, 20, 24, 25, 26, 16]).

It is worth noting that previous attention of MLAA reconstruction was focused on PET attenuation correction. In this paper, we exploit the MLAA reconstruction distinctly. We propose to combine the GCT image $\bm{\mu}$ with the X-ray CT image to obtain a dual-energy CT image pair to enable multi-material decomposition.

For each image pixel $j$, the high-energy GCT attenuation value $\mu_{j}$ and the low-energy X-ray CT attenuation value $x_{j}$ form a pair of dual-energy measurements $\bm{u}_{j}$. The tissue compositions are then described by a set of material bases, for example, air (A), soft tissue (S) or equivalently water, and bone (B):

$$\bm{u}_{j}\triangleq\left(\!\!\begin{array}[]{c}x_{j}\\ \mu_{j}\end{array}\!\!\right)=\underbrace{\left(\!\!\begin{array}[]{c c c c}x_{A}&x_{S}&x_{B}\\ \mu_{A}&\mu_{S}&\mu_{B}\end{array}\!\!\right)}_{\bm{U}}\underbrace{\left(\!\!\begin{array}[]{c}\rho_{j,A}\\ \rho_{j,S}\\ \rho_{j,B}\end{array}\!\!\right)}_{\bm{\rho}_{j}}$$

(5)

where the coefficients $\rho_{j,k}$ with $k=\{A,S,B\}$ are the fraction of each basis material in pixel $j$ and subject to

$$\sum_{k}\rho_{j,k}=1.$$

(6)

The material basis matrix $\bm{U}$ consists of the linear attenuation coefficients of each basis material measured at the low and high energies. The estimates of $\bm{\rho}_{j}$ are obtained using the following least-square optimization for each image pixel,

$$\hat{\bm{\rho}}_{j}=\arg\min_{\bm{\rho}_{j}}||\bm{u}_{j}-\bm{U}\bm{\rho}_{j}||^{2}.$$

(7)

The GCT by standard MLAA reconstruction is commonly noisy. To suppress the noise, we propose to utilize available X-ray CT in PET/CT as an image prior. As illustrated in figure 2, the higher contrast and potentially much better image quality provided by a X-ray CT image can be beneficial to guide the reconstruction of a $\gamma$-ray CT image. In this work, we apply the kernel method which was originally developed for dynamic PET reconstruction (e.g., [28]) and dual-modality imaging such as PET/MR (e.g., [29, 30]). Here we extend the kernel method to exploit low-energy X-ray CT image prior for reconstruction of 511 keV GCT image from PET emission data.

With the X-ray CT image $\bm{x}$, we can extract a “data point” $\bm{f}_{j}$ for each pixel $j$ from the image, for example, using the image patch centered at $j$ (Fig. 2). A transformed feature space can be defined by a nonlinear mapping function $\bm{\phi}$, which transforms the low-dimensional space $\{\bm{f}_{j}\}$ to a very-high dimensional space $\{\bm{\phi}(\bm{f}_{j})\}$. In the high-dimensional feature space, the intensity of the GCT in pixel $j$ can be described as a linear function, $\mu_{j}=\bm{w}^{T}\bm{\phi}(\bm{f}_{j}),$ where $\bm{w}$ denotes the coefficient vector. Because $\bm{w}$ also sits in the feature space, i.e., $\bm{w}=\sum_{l}\alpha_{l}\bm{\phi}(\bm{f}_{l})$, we then have the following equivalent kernel representation for $\mu_{j}$,

$$\mu_{j}=\sum_{l}\alpha_{l}\underbrace{\bm{\phi}(\bm{f}_{j})^{T}\bm{\phi}(\bm{f}_{l})}_{\kappa(\bm{f}_{j},\bm{f}_{l})}$$

(8)

where the kernel function $\kappa(\cdot,\cdot)$ is defined as the inner product of the two feature vectors $\bm{\phi}(\bm{f}_{j})$ and $\bm{\phi}(\bm{f}_{l})$. The form of the kernel function can be directly defined without knowing the specific form of $\bm{\phi}$. For example, the radial Gaussian kernel is

$$\kappa(\bm{f}_{j},\bm{f}_{k})=\exp\left(-||\bm{f}_{j}-\bm{f}_{k}||^{2}/2\sigma^{2}\right)$$

(9)

which corresponds to a $\bm{\phi}$ of infinite dimension. $\sigma$ is a hyper-parameter.

The matrix-vector form of the kernel representation for the GCT image is

$$\bm{\mu}=\bm{K}\bm{\alpha}.$$

(10)

where $\bm{K}$ is the kernel matrix built on the X-ray image $\bm{x}$ with its $(j,l)$th element equal to $\kappa(\bm{f}_{j},\bm{f}_{l})$. The unknown parameter vector $\bm{\alpha}$ denotes the corresponding kernel coefficient image. Although with a large matrix size, $\bm{K}$ can be built to be sparse to make a practical implementation [28].

Inserting the kernel representation in Eq. (10) into the original MLAA formulation leads to a kernelized optimization problem as follows,

$$\hat{\bm{\lambda}},\hat{\bm{\alpha}}=\arg\max_{\bm{\lambda}\geq\bm{0},\bm{\alpha}\geq\bm{0}}L\big{(}\bm{y}|\bm{\lambda},\bm{K}\bm{\alpha}\big{)}.$$

(11)

Once $\hat{\bm{\alpha}}$ is obtained, the final estimate of the GCT image is obtained by

$$\hat{\bm{\mu}}=\bm{K}\hat{\bm{\alpha}}.$$

(12)

To solve the optimization problem, we use the same alternating optimization strategy as used in [9]. Each iteration of the algorithm consists of two separate $\bm{\lambda}$-step and $\bm{\alpha}$-step,

	$\displaystyle\hat{\bm{\lambda}}$	$\displaystyle=$	$\displaystyle\arg\max_{\bm{\lambda}\geq\bm{0}}L\big{(}\bm{y}|\bm{\lambda},\bm{K}\hat{\bm{\alpha}}\big{)},$		(13)
	$\displaystyle\hat{\bm{\alpha}}$	$\displaystyle=$	$\displaystyle\arg\max_{\bm{\alpha}\geq\bm{0}}L\big{(}\bm{y}|\hat{\bm{\lambda}},\bm{K}\bm{\alpha}\big{)}.$		(14)

We simulated the GE Discovery 690 PET/CT scanner in 2D. The TOF timing resolution of this PET scanner was about 550 ps. The simulation was conducted using the XCAT phantom. The true PET activity image and 511 keV attenuation image are shown in Fig. 3(a) and (b), respectively. The images were first forward projected to generate noise-free PET sinogram of 11 TOF bins. A 40% uniform background was included to simulate random and scattered events. Poisson noise was then generated using 5 million expected events, unless specified otherwise. The x-ray CT image at a low-energy 80 keV is shown in Fig. 3(c). The data were reconstructed into images of $180\times 180$ with a pixel size of $3.9\times 3.9$ mm².

Three reconstruction algorithms were compared in this study: (1) the standard MLAA algorithm, (2) proposed kernel MLAA, and (3) post-reconstruction kernel smoothing using the same kernel matrix $\bm{K}$. The third algorithm is also equivalent to nonlocal means denoising [28]. Using the $3\times 3$ image patches extracted from the X-ray CT image $\bm{x}$, the kernel matrix was built using 50 nearest neighbors in a way similar to [28]. Each of the standard MLAA and kernel MLAA algorithms was run for 3000 iterations. Within each iteration, one inner iteration was used for the PET activity $\bm{\lambda}$ estimation step and five inner iterations were used for the attenuation $\bm{\mu}$ estimation step.

For each reconstruction algorithm, two different initial image estimates were used for the GCT reconstruction. One is the uniform initial with $\mu_{j}^{\mathrm{init}}=0.1~{}\textrm{cm}^{-1}$ and the other is the 511 keV attenuation map converted from the X-ray CT image using a bilinear scaling.

Our main interest is in dual-energy CT imaging. Hence we did not specifically evaluate the algorithms for PET activity reconstruction but focused on the evaluation for the CT performance.

The quality of GCT was first assessed using the image mean squared error (MSE) defined by

$$\mathrm{MSE}(\hat{\bm{\mu}})=10\log_{10}\frac{||\hat{\bm{\mu}}-\bm{\mu}^{\mathrm{true}}||^{2}}{||\bm{\mu}^{\mathrm{true}}||^{2}}\quad(\mathrm{dB}),$$

(29)

where $\hat{\bm{\mu}}$ is an image estimate of GCT obtained with one of the MLAA reconstruction methods and $\bm{\mu}^{\mathrm{true}}$ denotes the ground truth GCT image.

For evaluating quantification, we also calculated the ensemble bias and standard deviation (SD) of the mean intensity in regions of interest (ROI) by

	$\displaystyle\mathrm{Bias}$	$\displaystyle=$	$\displaystyle\frac{1}{c^{\mathrm{true}}}|\bar{c}-c^{\mathrm{true}}|,$		(30)
	$\displaystyle\mathrm{SD}$	$\displaystyle=$	$\displaystyle\frac{1}{c^{\mathrm{true}}}\sqrt{\frac{1}{N_{r}-1}\sum_{i=1}^{N_{r}}|c_{i}-\bar{c}|^{2}},$		(31)

where $c^{\mathrm{true}}$ is the noise-free regional intensity and $\bar{c}=\frac{1}{N_{r}}\sum_{i=1}^{N_{r}}c_{i}$ denotes the mean of $N_{r}$ realizations. $N_{r}=5$ in this study.

In addition to the comparison for GCT image quality, different reconstruction algorithms were further compared for dual-energy CT multi-material decomposition as formulated in Section 2. Image MSE, ROI bias and SD were calculated for each of the material basis fraction images.

Fig. 4 shows the reconstructed GCT attenuation images at 511 keV from the noisy PET emission data using the standard MLAA and proposed kernel MLAA algorithms with 400 iterations. Both the results of using the uniform initial and CT initial are shown. It is not surprising that the CT initial provided better image quality because the initial estimate is closer to the ground truth. For both initials, the kernel MLAA achieved much better results with lower MSE than the standard MLAA reconstruction.

Fig. 5 shows the resulting MSE as a function of iteration number in different reconstruction algorithms. The post-reconstruction denoising with kernel smoothing (KS) is also included in the comparison. For all the three reconstruction algorithms, the image initials made a large difference at early iterations but not at late iterations where the image reconstructions start to converge despite the initial starting point. In all the three reconstruction approaches, the CT initial also allowed an earlier iteration stopping to get each own best MSE than the uniform initial. This is useful as less number of iteration leads to accelerated speed. While post-reconstruction denoising improved the MLAA result, the kernel MLAA achieved a larger improvement on image quality with lower MSE. Because the CT initial demonstrated better performance than the uniform initial, hereafter we mainly present further comparisons based on the CT initial.

The effect of count level is shown in Fig. 6. In addition to the 5 million count level, two additional count levels (1 million and 10 million) were also included in the study. The number of iterations was fixed at 400 for each reconstruction. With increased count level, image quality by different algorithms were all increased. The kernel MLAA remains superior over the MLAA (with or without post-reconstruction smoothing) at different count levels.

Fig. 7 shows the results of ensemble bias versus SD for GCT ROI quantification in a liver region and a spine bone region. The count level was 5 million events. The iteration number varies from 0 to 3000 with a step of 100 iterations. As iteration number increases, the bias of ROI quantification is reduced while the SD is increased. After a certain number of iterations, the increasing noise may become dominant, which in turn induces higher bias. The post-reconstruction kernel smoothing approach outperformed the standard MLAA approach in a homogeneous region such as the liver but may oversmooth small targets such as the bone structures. The kernel MLAA achieved the best performance for both ROIs. At a fixed bias level, the kernel MLAA has lower SD than the other two approaches.

The results of applying multi-material decomposition (MMD) to the combined X-ray CT and GCT data are given in Figure 8. The MLAA and kernel MLAA reconstructions were run for 400 iterations. The ground truth of the three basis fractional images (air, soft tissue, bone) was generated using the noise-free pair of low-energy x-ray CT image and the 511 keV GCT image. The images by MLAA contain substantial noise but the regions of air and bone were still differentiated from the soft-tissue basis. Compared to MLAA, the kernel MLAA reconstruction led to a dramatic noise reduction in all the three basis images with increased image MSE.

The MSE of each basis fractional image is further plotted as a function of iteration number in Fig. 8(d) in which the post-reconstruction kernel smoothing approach was also included for comparison. The kernel MLAA approach demonstrated a significant MSE improvement over the conventional MLAA approach with or without post-reconstruction smoothing across all iterations.

To demonstrate the performance of different reconstruction algorithms for ROI quantification on MMD images, Fig. 10 shows the bias versus SD trade-off plot for ROI quantification on the bone fractional image using the spine ROI as shown in Fig. 7(a). Due to over-smoothing, the post-reconstruction denoising approach had lower SD but higher bias, resulting in an even worse trade-off than the MLAA without denoising. In comparison, the kernel MLAA reconstruction achieved a consistently better trade-off than the other two approaches.