Single-Subject Deep-Learning Image Reconstruction with a Neural Optimization Transfer Algorithm for PET-enabled Dual-Energy CT Imaging

Commonly the PET measurement $\bm{y}$ is well modeled as independent Poisson random variables which follows the log-likelihood function:

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

where $\bm{G}_{m}$ is the PET detection probability matrix for time-of-flight bin $m$. $\bm{r}_{m}$ accounts for the expectation of random and scattered events. $\bm{n}_{m}(\bm{\mu})$ is the normalization factor with the $i$th element being

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

The MLAA reconstruction algorithm [6] simultaneously estimates the attenuation image $\bm{\mu}$ and the activity image $\bm{\lambda}$ from the PET projection data $\bm{y}$ by maximizing the Poisson log-likelihood,

An iterative interleaved updating strategy is commonly used to seek the solution. At each iteration of the algorithm, $\bm{\lambda}$ is first obtained based on the attenuation image $\bm{\mu}^{n}$ from the previous iteration $n$:

which can be updated by the maximum-likelihood expectation maximization (MLEM) algorithm [29] with one subiteration,

where $\bm{p}^{n}$ denotes the updated sensitivity image from the iteration $n$, $\bm{p}^{n}=\sum_{m}\bm{G}_{m}^{T}\bm{n}_{m}(\bm{\mu}^{n}).$

$\bm{\mu}$ is then updated with the estimated $\bm{\lambda}$ following the maximum-likelihood transmission reconstruction formulation,

where

with $\bm{l}=\bm{A}\bm{\mu}$ and

The sub-optimization problem Eq. (7) can be solved using the separable paraboloidal surrogate algorithm [30].

Note that conventional applications of MLAA mainly focused on improving PET attenuation correction (e.g., [31, 32, 8, 33, 34, 35, 36]). A new application of MLAA demonstrated the potential of gCT reconstruction for improving the estimation of proton stopping-power in proton therapy [37]. Differently in our PET-enabled DECT method, the gCT image $\bm{\mu}$ is combined with X-ray CT image $\bm{x}$ to form a DECT image pair for multi-material decomposition [5].

The gCT image estimate by the MLAA method [6] is commonly noisy due to the limited counting statistics of PET emission data. To suppress noise, the kernel MLAA or KAA integrates the X-ray CT prior image into the PET forward model by describing the gCT image intensity $\mu_{j}$ at pixel $j$ using a kernel representation [5, 38],

where $\kappa(\cdot,\cdot)$ is the kernel function (e.g., radial Gaussian) with $\bm{f}_{j}$ and $\bm{f}_{l}$ denoting the feature vectors of pixel $j$ and $l$ that are extracted from the X-ray CT image $\bm{x}$. $\mathcal{N}_{j}$ defines the neighborhood of pixel $j$, for example, selected by a k-nearest neighbor algorithm. $\alpha_{l}$ denotes the corresponding kernel coefficient of each neighboring pixel $l$ in $\mathcal{N}_{j}$.

The equivalent matrix-vector form for the gCT image representation is

where $\bm{K}$ is the kernel matrix and $\bm{\alpha}$ is the corresponding kernel coefficient image. Substituting Eq. (11) into the MLAA formulation in Eq. (4) gives the following KAA optimization formulation,

Once $\hat{\bm{\alpha}}$ is determined, the gCT image is obtained as $\hat{\bm{\mu}}=\bm{K}\hat{\bm{\alpha}}$. Note that the conventional MLAA can be considered as a special case of the KAA with $\bm{K}$ equal to an identity matrix.

The KAA problem is also solved using an interleaving optimization strategy between the activity image $\bm{\lambda}$ update and the kernel coefficient image $\bm{\alpha}$ update [5]. In each iteration of KAA, $\bm{\lambda}^{n+1}$ is first obtained using the MLEM updating formula Eq. (6) and then $\bm{\alpha}^{n+1}$ is obtained using the following kernelized transmission reconstruction (KTR) optimization,

as illustrated in Fig. 2a.

For each image pixel $j$, the gCT attenuation value $\mu_{j}$ and X-ray CT attenuation value $x_{j}$ jointly form a pair of dual-energy measurements $\bm{u}_{j}\triangleq[x_{j},\mu_{j}]^{T}$, which can be modeled by a set of material bases, such as air (A), soft tissue (S) or equivalently water, and bone (B):

subject to $\sum_{k}\rho_{j,k}=1.$ The coefficients $\rho_{j,k}$ with $k={A,S,B}$ are the fraction of each basis material in pixel $j$. The material basis matrix $\bm{U}$ consists of the linear attenuation coefficients of each basis material measured at the low and high energies. Finally, $\bm{\rho}_{j}$ is estimated using the following least-square optimization for each image pixel,

While demonstrating a substantially better performance than MLAA (e.g., in [5]), KAA may still suffer from noise or artifacts in low-count cases. In this work, we exploit neural networks as a conditional deep coefficient prior for improving KAA for gCT image reconstruction.

The kernel coefficient image $\bm{\alpha}$ for $\bm{\mu}$ in Eq. (11) is described as a function of the conditional neural networks,

where $\bm{z}$ is the available image prior from the corresponding X-ray CT in this work and $\bm{\psi}$ denotes the neural network mapping from the known input image $\bm{z}$ to the $\bm{\alpha}$ image with $\bm{\theta}$ the parameters of the neural network.

The gCT image is then modeled using the following kernel representation with the conditional deep coefficient prior,

Fig. 3 shows a graphical illustration of the proposed model for representing a gCT image, of which the last layer is a linear kernel representation with pre-determined weights $\{\kappa_{j,l}\}$ that are also calculated from $\bm{z}$.

By substituting the gCT representation Eq. (17) into the forward model in Eq. 2, we have the following forward projection model that is parameterized by the neural network parameters $\bm{\theta}$:

The model is equivalent to the standard KAA image model in Eq. (2) if the neural network $\bm{\psi}$ is an identity mapping. It is also equal to a conditional DIP (CDIP) model directly in the gCT image domain if the kernel matrix $\bm{K}$ is an identity matrix, which leads to $\bm{\mu}=\bm{\psi}(\bm{\theta}|\bm{z})$.

By combining the forward projection model with the MLAA formulation Eq. (4), we have the proposed neural KAA optimization formulation,

Similar to the optimization approach for MLAA and KAA, an interleaving updating strategy can be used here to estimate $\bm{\lambda}$ and $\bm{\theta}$ iteratively,

Once $\bm{\theta}$ is estimated, the gCT image is obtained by

Compared to the previous KAA approach, the neural KAA approach combines the kernel representation with a neural network-based deep coefficient prior, which introduces an implicit regularization for KAA to improve the gCT image estimate $\hat{\bm{\mu}}$ through Eq. (21).

In each iteration of the neural KAA, the $\bm{\lambda}$ estimation step in Eq. (20) can be directly implemented by using the MLEM algorithm Eq. (6). The $\bm{\theta}$-estimation step Eq. (21) follows a KTR formulation but with using neural networks as a conditional deep coefficient prior. For simplicity, we rewrite this neural KTR problem in Eq. (21) as

where $J_{n}(\bm{\theta})$ is the transmission likelihood function at iteration $n$,

with $h_{i,m}$ defined in Eq. (8) but the line integral $l$ following a kernelized neural-network model,

The optimization problem Eq. (23) is challenging to solve because the unknown neural network parameters $\bm{\theta}$ is nonlinearly involved in the projection domain for transmission imaging due to Eq. (25), resulting in a complex optimization problem.

One commonly used solution would be a type of gradient descent algorithm that uses the chain rule to calculate the gradient of $J_{n}(\bm{\theta})$ with respect to $\bm{\theta}$ [19, 17, 18, 25],

which is then fed into an existing deep learning package to estimate the network parameters from the PET projection data $\bm{y}$. However, such an approach ties each calculation of $\frac{\partial\bm{\psi}}{\partial\bm{\theta}}$, which relates to the neural network component, with a calculation of $\frac{\partial J_{n}(\bm{\theta})}{\partial\bm{\psi}}$, which relates to the PET tomographic reconstruction. While the former operation is usually efficient by using a deep learning library, the latter operation requires both forward and back projections of PET data and is computationally intensive due to the large size of the transmission system matrix $\bm{A}$. As a result, the whole algorithm can be slow due to the natural need for many iterations for training a neural network.

Another solution would be the ADMM algorithm as used for the DIP reconstruction [13, 20, 21]. This type of algorithms may separate the neural network learning step (that uses $\frac{\partial\bm{\psi}}{\partial\bm{\theta}}$) from the tomographic reconstruction step (that uses $\frac{\partial J_{n}(\bm{\theta})}{\partial\bm{\psi}}$). However, a major weakness of ADMM is that it is difficult to tune the induced hyper-parameters in practice as demonstrated by many studies [22, 23].

In this work, we develop an approach to decoupling the reconstruction step and neural network learning step without adding extra hyper-parameters towards a more practical and efficient implementation.

The neural network reconstruction problem shares the same complication as the tomographic reconstruction of nonlinear kinetic parameters in dynamic PET [40] which, however, can be solved by a quadratic surrogate-based optimization transfer algorithm [39]. Thus we apply the principle of optimization transfer [24] to convert the original difficult problem for $J_{n}(\bm{\theta})$ into a quadratic surrogate optimization problem, as illustrated in Fig. 4. The quadratic surrogate function $\phi_{J}(\bm{\theta}|\bm{\theta}^{n})$ minorizes the transmission likelihood function $J_{n}(\bm{\theta})$,

Then the maximization of $J_{n}(\bm{\theta})$ is transferred into maximizing the surrogate function,

The quadratic surrogate $\phi_{J}(\bm{\theta}|\bm{\theta}^{n})$ is designed to be easy to solve. The new update $\bm{\theta}^{n+1}$ is guaranteed to monotonically increase the original likelihood $J_{n}(\bm{\theta})$, i.e.,

as also demonstrated by Fig. 4.

Note that the concept of optimization transfer has been also applied for the KAA reconstruction in which the gCT projection is linear with respect to the unknown attenuation image $\bm{\mu}$ [5]. However, in the neural KAA, the unknown becomes the neural network parameters $\bm{\theta}$ which is non-linearly involved in the gCT projection domain. Our derivations will show how the quadratic surrogate functions are built to decouple the linear tomographic reconstruction step and the nonlinear neural network learning step from each other, as also shown in Fig. 2.

One difficulty with dealing the objective function $J_{n}(\bm{\theta})$ directly is the Poisson log-likelihood function follows a non-quadratic form. Following a derivation similar to [41, 30], a paraboloidal surrogate function can be constructed for $J_{n}(\bm{\theta})$,

where $l_{i}^{n}$ is $l_{i}$ calculated with $\bm{\theta}^{n}$ and $\Delta l_{i}=l_{i}-l_{i}^{n}$. $\eta_{i,m}(l)$ is chosen by design as the optimal curvature of the Poisson log-likelihood [41],

where $\dot{h}_{i,m}(l)$ and $\ddot{h}_{i,m}(l)$ are the first and second derivatives of $h_{i,m}(l)$, respectively [41].

With several algebraic operations similar to the way in [39], we then can derive the following equivalent quadratic form for $S(\bm{\theta}|\bm{\theta}^{n})$ in Eq. (31),

where $C_{S}^{n}$ is the remainder that is independent of the unknown parameter $\bm{\theta}$. $\hat{\bm{l}}^{n+1}$ is an intermediate gCT projection data [5],

and $\hat{\bm{\eta}}^{n}$ is an intermediate weight also in the projection domain,

Based on this quadratic surrogate $S(\bm{\theta}|\bm{\theta}^{n})$, the maximization of $J_{n}(\bm{\theta})$ for the neural KTR is transferred to an equivalent weighted least-square reconstruction problem,

This quadratic form is simpler than the original non-quadratic likelihood function $J_{n}(\bm{\theta})$. However, the nonlinear neural network model $\bm{\psi}(\bm{\theta}|\bm{z})$ is still coupled in the projection domain.

By considering $\bm{A}\bm{K}$ as a single matrix and using the convexity of Eq. (33) to build a separable quadratic surrogate [42], we can construct the following surrogate $s(\bm{\theta}|\bm{\theta}^{n})$ for $S(\bm{\theta}|\bm{\theta}^{n})$,

where $\Delta\bm{\psi}(\bm{\theta}|\bm{z})\triangleq\bm{\psi}(\bm{\theta}|\bm{z})-\bm{\psi}(\bm{\theta}^{n}|\bm{z})$. $\bm{g}^{n}$ is the gradient of $S(\bm{\theta}|\bm{\theta}^{n})$ with respect to $\bm{\psi}$ at iteration $n$,

and $\bm{\omega}^{n}$ is an intermediate weight image,

where $\mathbf{1}$ is the all-one vector. Note that $\bm{g}^{n}$ is also equal to the gradient of $J_{n}(\bm{\theta})$ with respect to $\bm{\psi}$ at iteration $n$.

After matching a quadratic function, we can have the following equivalent form for $s(\bm{\theta}|\bm{\theta}^{n})$,

where $\hat{\bm{\alpha}}^{n+1}$ is an intermediate kernel coefficient image calculated from

with $\bm{\alpha}^{n}\triangleq\bm{\psi}(\bm{\theta}^{n}|\bm{z})$. $C_{s}^{n}$ denotes the corresponding remainder term that is independent of $\bm{\theta}$. Note that $\hat{\bm{\alpha}}^{n+1}$ is equivalent to one iteration of the KTR reconstruction algorithm in [5].

Thus, the weighted least-square reconstruction problem Eq. (36) defined in the projection domain is transferred to the following neural-network learning problem that is fully defined in the image domain,

In particular, the cost function here follows a weighted $l_{2}$-norm loss that can be easily used with an existing deep learning package (e.g., PyTorch).

It is straightforward to prove that the surrogate function $s(\bm{\theta}|\bm{\theta}^{n})$ minorizes the original likelihood function $J_{n}(\bm{\theta})$,

Note that the optimization transfer algorithm in [5] for standard KAA (Fig. 2a) is only a special case of the derived algorithm for the neural KAA (Fig. 2b) here without having the least-square neural network learning step. Therefore we call the new algorithm the neural optimization transfer algorithm to highlight its estimation of the neural network model parameters. The associated least-square form in Eq. (42) for network training is not arbitrarily defined but analytically derived from the theory of optimization transfer for a non-linear model. Furthermore, no any new hyper-parameters are introduced from the optimization transfer process, which represents an important advantage as compared to a possible ADMM algorithm.

A pseudo-code of the proposed neural optimization transfer algorithm for gCT reconstruction is summarized in Algorithm 1. The algorithm consists of three separate steps in each iteration:

1. 1.

   PET activity image reconstruction: Obtain the PET activity image update $\bm{\lambda}^{n+1}$ using one iteration of MLEM algorithm in Eq. (6).

2. 2.

   Kernelized gCT image reconstruction: Obtain an intermediate kernel coefficient image update $\hat{\bm{\alpha}}^{n+1}$ using one iteration of the KTR algorithm in Eq. (41);

3. 3.

   Least-square neural-network learning for gCT image approximation: Update the network model parameters $\bm{\theta}^{n+1}$ using Eq. (42) to approximate $\hat{\bm{\alpha}}^{n+1}$;

Step 1 and step 2 are updated analytically. The tomographic reconstruction of projection data are only involved in these two steps. Step 3 can be implemented efficiently using existing deep learning packages without involving any projection data directly. Thus, the neural network learning step is decoupled from the image reconstruction steps and is easy to implement in practice. The initialized weights of the neural network learning at the reconstruction iteration $n+1$ are inherited from the previous iteration $n$ rather than a random initialization. Note that random initialization was also tested but demonstrated a less stable performance. The algorithm is guaranteed to increase the likelihood function monotonically due to the nature of optimization transfer [24].

Note that many other imaging modalities such as X-ray CT [17, 20], MRI [25, 21], and optical tomography [26, 27] may often use the same weighted least-square reconstruction formula in Eq. (36) (but with constant $\hat{\bm{l}}^{n+1}$ and $\hat{\bm{\eta}}^{n}$). Thus the neural optimization transfer algorithm in Algorithm 1 may be also applicable to these reconstruction problems, though without the need for Step 1 which is specific for the neural KAA for gCT reconstruction.

Any neural-network model that is suitable for image representation would be compatible with the proposed algorithm. As an example, here we use a popular residual U-Net architecture to represent the kernel coefficient image, as illustrated in Fig. 5, following its previous successful use for PET activity image reconstruction [13, 43]. The model comprises the following key components: (1) A convolution layer (kernel size: 3) + batch normalization (BN) + leaky rectified linear unit (LReLU) for feature extraction; (2) The strided convolution layer (strided size: 2, kernel size: 3) for down-sampling; (3) Bilinear/trilinear interpolation for up-sampling; (4) The additive operation used for skip connection between encoder and decoder paths; (5) A ReLU operation prior to the output layer used for generating the non-negative kernel coefficient image. The number of feature maps for each layer is also included in Fig. 5.

In this work, the input image of the neural network was set to an X-ray CT image, that can be explained as the conditional deep image prior (CDIP) [44]. We used the Adam algorithm with a learning rate of $10^{-3}$ and 150 sub-iterations for the least-square neural network learning, which were the same setting as previous work [43]. The neural-network learning step was implemented with PyTorch on a PC with an Intel i9-9920X CPU with 64GB RAM and a NVIDIA GeForce RTX 2080Ti GPU.

We first conducted a two-dimensional computer simulation study following the GE Discovery 690 PET/CT scanner geometry. This PET scanner has a time-of-flight (TOF) timing resolution of approximately 550 ps. The simulation was conducted using one chest slice of the XCAT phantom [45]. Fig. 6a and Fig. 6b show the simulated ground truth of PET activity image and 511 keV gCT attenuation image, respectively. The low-energy X-ray CT image was simulated from XCAT at 80 keV and is shown in Fig. 6c. The activity and gCT images were first forward projected to generate noise-free emission sinograms with the size of 281 (radial bins) $\times$ 288 (angular bins) $\times$ 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. The projection data was reconstructed into PET activity and gCT images of 180$\times$180 with a pixel size of 3.9$\times$3.9 mm². Ten noisy realizations were simulated and reconstructed for comparing reconstruction methods.

In this work, four types of reconstruction methods were compared, including (1) the standard MLAA [6], (2) existing KAA [5], (3) a CDIP reconstruction method, which is equal to the proposed neural KAA with $\bm{K}=\bm{I}$, and (4) proposed neural KAA that uses the neural optimization transfer algorithm.

For constructing kernels, the feature vector $\bm{f}_{j}$ was chosen as the pixel intensities of X-ray CT image in a $3\times 3$ image patch centered at pixel $j$. The radial Gaussian kernel function, $\kappa(\bm{f}_{j},\bm{f}_{l})=\exp(-||\bm{f}_{j}-\bm{f}_{l}||^{2}/2\sigma^{2})$, was used to build the kernel matrix $\bm{K}$ using $\sigma=1$ and $k$NN search with $k=50$, in the same way as [5].

The initial estimate of the PET activity image was set to a uniform image. Following [5], we used the X-ray CT-converted 511 keV attenuation map as the initial estimate of gCT image for accelerated convergence. All reconstructions were run for 3000 iterations for investigating the convergence behaviors of different algorithms.

As the focus of this work is on gCT for DECT, the evlauation of PET activity image quality is not concerned. Different reconstruction algorithms were first compared for gCT image quality using the mean squared error (MSE),

where $\hat{\bm{\mu}}$ represents the reconstructed gCT image by each method and $\bm{\mu}^{\rm{true}}$ denotes the ground truth. The error image, defined as $\hat{\bm{\mu}}-\bm{\mu}^{\rm{true}}$, was used for highlighting the differences in key parts. The ensemble bias and standard deviation (SD) of the mean intensity in a regions of interest (ROI) were also calculated to evaluate ROI quantification in a liver region and a bone region (Fig. 6(d)),

where $c^{\rm{true}}$ is the true average intensity in a ROI and $\overline{c}=\frac{1}{N_{r}}\sum_{i=1}^{N_{r}}c_{i}$ denotes the mean of $N_{r}$ realizations ($N_{r}=10$).

Different reconstruction algorithms were further compared for DECT multi-material decomposition. Similarly, the image MSE, error image, and ROI-based bias and SD were calculated for each of the material basis fractional images.

Fig. 7 shows the plots of Poisson log-likelihood as a function of iteration number for one simulated data for the proposed neural optimization transfer and gradient descent algorithms. Here the implementation of gradient descent algorithm followed [17] with an optimized step size using grid search. While the gradient descent algorithm not surprisingly demonstrates an oscillating behavior because of the non-convex nature of the original likelihood function (Eq. 23), the proposed neural optimization transfer algorithm shows a monotonic increase in the likelihood function as the iteration increases, which is guaranteed by the theory of optimization transfer.

Fig. 24 shows the true and reconstructed 511-keV gCT images using different algorithms with 400 iterations, as well as corresponding error images. The image MSE results were included for quantitative comparison. The selection of 400 iterations was a balance for different methods (KAA, CDIP and the proposed method) to show good image quality according to the MSE performance. The MLAA reconstruction was noisy. The KAA method substantially improved the gCT image quality, but with a lower contrast in the bone region as compared to the ground truth (as evidenced in its error image). The CDIP method had a slightly better MSE than KAA, but it induced artifacts, which was in turn propagated into the material decomposition results as shown later in Fig. 25. In comparison, the proposed neural KAA demonstrated the least level of noise with good visual quality, and achieved the lowest MSE among different algorithms.

Fig. 9a shows the MSE plots as a function of iteration number for different algorithms. The iteration number varied from 0 to 3000 with a step of 100 iterations and error bars were calculated over the 10 noisy realizations. Compared to KAA, the proposed neural KAA showed a lower MSE among different methods. The curve also shows that similar to other algorithms, early stopping of the iterations is beneficial for the neural KAA to obtain good image quality while keep computational efficiency.

Fig. 9b and Fig. 9c further show the comparison of ensemble bias versus SD for gCT quantification in a liver region and a bone region. The curves were obtained by varying the iteration number from 300 to 3000 iterations with an interval of 100 iterations. As iteration number increases, the bias of ROI quantification is reduced while the SD increases. At a comparable bias level, the proposed neural KAA had a lower noise SD than the other three methods.

We also conducted a comparison of different reconstruction methods for multi-material decomposition (MMD). Fig. 25 shows the fractional basis images of bone and soft tissue obtained from MMD of the PET-enabled DECT images with 400 iterations, as well as their corresponding error images. The ground truth of the soft tissue and bone bases was generated using the noise-free pair of low-energy X-ray CT image ($\bm{x}$ in Eq. (14)) and the 511 keV gCT image. The conventional KAA method outperformed MLAA but were still with noise artifacts, as pointed by arrows. Interestingly, even though the CDIP reconstruction had a better MSE for gCT images than KAA, the benefit did not propagate into the MMD basis images. The proposed neural KAA achieved less noise and artifacts than KAA thanks to the CDIP-based regularization on the kernel coefficient image $\bm{\alpha}$.

Fig. 11 further shows image MSE as a function of iteration number for each basis fractional image. Again, the proposed neural KAA demonstrated the lower minimum MSE result in each basis image among different reconstruction methods.

Fig. 12 shows the quantitative comparisons of ensemble bias versus SD for ROI quantification on the soft tissue (a) and bone (b) fractional images by varying the reconstruction iteration number. Similar to the results of gCT ROI quantification, the neural KAA achieved the lowest noise level at a comparable bias level. It is noticeable that for the bone ROI quantification, both the standard KAA and neural KAA showed a bias when compared to MLAA. The bias was propagated from the gCT reconstruction (as shown in Fig. 9c).

Our experiments indicated that neural network learning is stable when the learning rate in the Adam optimizer ranges from $10^{-4}$ to $10^{-2}$. A larger rate may make the learning difficult to converge, while a smaller rate may reduce the convergence rate. The sub-iteration number used in the least-square neural network learning may also have influenced the results. Fig. 13 shows the effect of this sub-iteration number of neural network learning and total outer-iteration number of the neural KAA on the gCT image MSE. The results suggest the algorithm becomes relatively stable after 150 sub-iterations under different total outer-iterations of the neural KAA.

We have further evaluated different reconstruction methods using a real three-dimensional (3D) phantom scan on the uEXPLORER PET/CT scanner [46] at UC Davis. This phantom [47] was filled with water in the background and four inserts were filled with (1) lung tissue equivalent material, (2) water, and (3&4) salt water, as shown in Fig. 14. Other attenuation materials composed of fat tissue-equivalent materials and bovine rib bones were wrapped around the phantom. A ¹⁸F-FDG solution was uniformly filled in all five compartments. Two X-ray CT scans, one at 80 kVp and the other at 140 kVp, were acquired to provide the reference for MMD analysis. In this study, the reconstructions were performed on a truncated dataset of 2-min scan duration and 7 cm axial length to reduce computational time. The projection data were acquired with 27 TOF bins for a TOF resolution of 505ps [46]. The 3D sinogram for each TOF bin was of 533 radial bins, 420 angular bins, and 576 direct and oblique detector planes. The reconstructed gCT image size was $150\times 150\times 17$ with a voxel size of $4\times 4\times 4$ mm³. For MMD using a PET-enabled DECT image pair, the X-ray CT image was downsampled to match with the gCT voxel size using linear interpolation.

We used the 80 kVp CT image as the image prior to generate the kernel matrix $\bm{K}$ and the input of the neural network. The settings of the neural network in the CDIP method and proposed neural KAA were the same as in the simulation study except using a 3D version of U-Net to match with the data. All reconstructions were implemented using the CASToR package [48] as described in [47] and run for 400 iterations. The 80 kVp X-ray CT-converted attenuation map was used as the initial for gCT.

Fig. 15 shows the gCT images reconstructed using different algorithms with 400 iterations. Similar to the results shown in the simulation, the MLAA was extremely noisy. Both the standard KAA and CDIP methods reduced noise significantly but still contained significant artifacts. In contrast, the proposed neural KAA demonstrated the best visual quality.

Fig. 16 shows the bone fractional images from MMD and corresponding error images using different approaches. Compared to the reference from X-ray DECT, the image by MLAA demonstrated heavy noise. Both the KAA and CDIP methods suppressed the noise, but not without additional noise or artifacts. In comparison, the proposed neural KAA demonstrated the least artifacts and noise in the uniform regions and achieved the most similar bone fraction pattern with the reference image, as pointed by the arrows.

Fig. 17 further shows a quantitative comparison for ROI quantification of bone fraction. Here the ROI quantification is plotted versus the background noise SD measured in the water region by varying the iteration number from 40 to 400 with an interval of 20 iterations. As the iteration number increases, the estimated bone fraction becomes lower while the SD increases. Compared to the other three methods, the proposed neural KAA was the closest to the X-ray DECT reference (dashed line) and also achieved the lowest background noise level.

Based on the maximum-a-posteriori (MAP) reconstruction of PET emission data, we can update the attenuation image $\bm{\mu}$ following

where $\bm{R}(\bm{\mu})$ is a regularization term that imposes prior information on $\bm{\mu}$, controlled by the regularization parameter $\beta$. By using the FBS algorithm as the same as used in [3], the optimization of Eq. (47) is performed in the following steps

where Eq. (48) is a gradient descent step of $\bm{R}$ with the step size of $\gamma$, which can be solved by a deep-learning model, e.g., U-Net. Eq. (49) is an optimization problem regarding the log-likelihood function with $1/\gamma$ as the regularization hyperparameter. Different from [3] for PET activity image reconstruction, here our focus is for gCT image reconstruction and the separable paraboloidal surrogate (SPS) algorithm is thus used

where $\bm{\mu}_{\rm{SPS}}^{n+1}$ is updated by

where $\bm{g}^{n}$ and $\bm{\omega}^{n}$ are the gradient image and intermediate weight image respectively as given by same formulas to Eq. (38) and Eq. (39) with an identity matrix $\bm{K}$ and an identity mapping $\bm{\psi}$.

Finally, by setting the derivative of Eq. (50) to zero, we get a closed-form solution

We call this new algorithm a transmission-reconstruction or FBSTR to emphasize its nature for gCT attenuation image reconstruction.

The unrolled FBSTR algorithm can be represented as a deep iterative neural network, FBSTR-Net $\bm{\phi}_{FBSTR}$, with $N$ blocks, as illustrated in Fig. 21. For each block, we have three steps: (1) the regularization update based on the previous gCT image estimate using a U-Net model; (2) MLTR from projection domain using the SPS (Eq. (51)); (3) pixel-by-pixel image fusion by Eq. (52). The trainable parameters of U-Net model are shared across all blocks. We formulate the training process between the output of FBSTR-Net ($\bm{\mu}_{s}^{N}$) and the reference image ($\bm{\mu}_{s}^{\rm{Ref}}$) using a mean-squared-error (MSE) loss function

with $\bm{\mu}_{s}^{N}=\bm{\phi}_{FBSTR}(\bm{\theta}|\bm{y}_{s},\bm{\mu}^{0}_{s})$. $N_{s}$ is the number of training data. $\bm{\theta}$ includes the trainable parameters in the U-Net model and the hyperparameter $\gamma$. $\bm{y}_{s}$ is the projection data and $\bm{\mu}^{0}_{s}$ is the initial estimate.

Following the same simulation setting as described in Section V.A, we established a training dataset using ten other chest slices of the XCAT phantom, following a similar strategy as used in [2]. The examples of gCT images for training are shown in Fig. 22. For each case, we simulated ten noisy realizations so that total 100 training samples were included in our experiment. The ratio of the training set to the validation set is 2 to 8 [3]. The simulated data in Section V.A was considered as the testing data for method comparison. The same U-Net structure used for the neural KAA and CDIP was used in the FBSTR-Net instead of the residual learning unit used in [3] because of the U-Net’s better performance. The block number $N$ and learning rate was selected as 50 and 0.01 respectively according to the performance of validation set. The initialization, $\bm{\mu}^{0}$, was also the X-ray CT image-converted 511 keV attenuation map.

Fig. 23 shows the plots of training loss and validation loss with the change of training epoch number. Here the loss value (Eq. (53)) is converted into MSE measurement in dB. We thus chose the trained model at 150 epochs for testing.

Fig. 24 shows the true and reconstructed 511-keV gCT images using different algorithms as well as the corresponding error images. Here, we mainly focus on the FBSTR-Net results. The image MSE value was also included for quantitative comparison. It can be seen that the gCT image derived from FBSTR-Net had a slightly worse MSE than KAA. Its error image demonstrated relatively larger bias around bone regions. The ensemble bias and SD for gCT quantification in a liver region and a bone region using different reconstruction methods are presented in Table I. In addition to FBSTR-Net, the quantification results of 400 iterations of other methods were also included for comparison. A larger bone bias was observed in the FBSTR-Net result, which can be reflected from the error image in Fig. 24.

Fig. 25 shows the fractional basis images of bone and soft tissue obtained from MMD of the PET-enabled DECT images. The corresponding error images are also shown. The FBSTR-Net demonstrated a superior performance than CDIP and was comparable to KAA. While the CDIP had a better MSE than FBSTR-Net for gCT image comparison in Fig. 24, in the case of fractional images in Fig. 25, the inferior MSE results of CDIP were attributable to significant artifacts in the uniform regions. In comparison, the proposed neural KAA achieved the best result with the lowest MSE. Compared to KAA and FBSTR-Net, the improvement is due to the additional use of neural network as the deep coefficient prior ($\bm{\alpha}$). Similar to the Table I, the quantitative comparisons of ensemble bias and SD for ROI quantification on the soft tissue and bone fractional images are presented in Table II. While FBSTR-Net exhibited a bias in the soft tissue region comparable to that of KAA, it demonstrated a more pronounced bias in the bone region, which was propagated from the gCT reconstruction.

By and large, the FBSTR-Net, an example of model-based deep reconstruction methods, suggested a promising direction for gCT image reconstruction. Once trained, it can be applied directly to other data. However, it is challenging to apply the method to patient scans directly given no training database has been built. In comparison, the single-subject deep learning reconstruction method is directly applicable to a patient scan. In the future, we will apply the population-based deep learning approach for the real data when a training dataset of PET-enabled DECT imaging becomes available and investigate its generalization performance on unseen PET emission data.

The architecture of neural networks is crucial for training a model-based deep reconstruction method. In this preliminary study, we only compared the U-Net architecture and a general residual learning unit (no downsampling) [3]. A better performance may be achieved by using a more advanced neural-network architecture. When applying such model-based methods, it will be important to also elaborately investigate the model hyperparameters (e.g., block number, learning rate, optimizer) and training strategies. A comprehensive investigation of these aspects for FBSTR-Net is beyond the scope of this paper but would be worth exploring in the future.