2.1.

XCAT Phantom Definition

Voxelized phantoms were generated using the XCAT phantom program developed by Segars et al. at Duke University.²¹ Because many cancer indications relevant for treatment with proton therapy involve tumors in, or in close proximity to, the brain, spinal cord, or other constituents of the central nervous system, the head was chosen as the focus region in this work. The phantom anatomy used defined the skull and brain region, with skull and skin thicknesses set to 6.5 and 2.5 mm, respectively.²²

To model anatomical variation, six different XCAT phantoms were included in the data set (three male and three female). One male (“vmale50”) and one female phantom (“vfemale50”) were of standard 50th percentile geometries while the remaining phantoms (“male_pt77,” “male_pt148,” “female_pt86,” and “female_pt147”) introduced further variation. The male phantoms had BMIs of 22.71, 24.38, and $32.18\mathrm{kg}/{\mathrm{m}}^2$, and the female phantoms 18.21, 24.06, and $38.81\mathrm{kg}/{\mathrm{m}}^2$. However, because this work focused on the head, neither BMI or gender was expected to contribute significantly to the variation, and a more meaningful measure is the range of head sizes included. The maximal geometrical diameters of the phantoms were 203, 208, 209, 215, 223, and 224 mm, and maximal water equivalent thicknesses were 187, 193, 193, 200, 204, and 206 mm.²³ The anatomical variation of the XCAT phantoms is generated from the PeopleSize program, based on anthropometric dimensions of individuals from several countries.²⁴ From each phantom, 281 slices were selected, yielding a data set of 1686 examples.

The XCAT program includes the functionality of “activity” phantom generation, originally intended for nuclear medicine simulations. The output is a three-dimensional matrix where voxels of a delineated body structure are assigned the same number. This feature was used to assign each voxel a known tissue material, with compositions in accordance with Appendix A from ICRU Report 46.²⁵

To introduce substructure in the brain, the brain structures were classified as white or gray matter. If the structure contained a mix, it was assigned the whole brain material from Ref. 25. The white and gray matter materials were elementally defined according to densities and compositions listed in Ref. 26. The XCAT phantom included a background structure, comprising the space between the brain and the skull, as well as between the skull and the skin. The background was set to adipose, in accordance with pre-set values for the attenuation phantom generation by the XCAT program.

The phantom pixel width and slice width were set to 0.25 and 0.416 mm, respectively. This choice was motivated by the desire to keep the phantom resolution similar to or higher than the size of the reconstructed pixels. In addition, using a slice thickness equal to that of the reconstructed images facilitated the pairing of corresponding CT and ground truth SPR slices.

2.2.

Ground Truth SPR Maps

Based on the known material composition in each voxel, the XCAT phantom was converted to ground truth SPR maps. The relative electron density and mean ionization energy of the materials were calculated as

Three soft tissue materials (white matter, gray matter, and whole brain) and the skull were selected for evaluation of SPR accuracy. Three regions of interest (ROI) for each material were selected in one slice. RMSE, relative error, and relative standard deviation were calculated as percentages for each ROI, with formulas given by

Fig. 1

ROIs of the (a) white matter, (b) gray matter, (c) whole brain, and (d) skull materials. For all materials, ROIs 1, 2, and 3 are drawn as red, blue, and green, respectively.

Table 1

Number of pixels in the ROIs.

2.3.

Photon-Counting CT Simulations

Photon-counting CT imaging of the XCAT phantom was performed in CatSim²⁹ using a proprietary simulation model of a deep-silicon-based photon-counting CT prototype system developed by GE HealthCare.³⁰ An open-source version of CatSim, with more limited functionality, is available online.³¹

The simulations were defined by parameters corresponding to a realistic CT protocol for a head scan. The tube current, rotation time, tube voltage, and number of view angles were set to 260 mA, 1 s, 120 kV, and 4000, respectively. We simulated an axial scan with 5 mm coverage, using an approximation of the GE LightSpeed large bowtie filter (“VCTlarge.txt”). For each view and detector pixel, a single projection line was simulated.

Physical phenomena that are known to impact image quality can be included in the simulations, such as electronic noise (manifested as Poisson-distributed counts in the lowest energy bin), quantum noise, pileup (photons arriving at the detector during the dead time following the registration of a previous photon, which consequently go unregistered and/or distort the detected energy spectrum), and crosstalk (charge sharing and reabsorption of Compton scattered photons). The simulations in this study included electronic and quantum noise, but excluded the effects of crosstalk and pileup, to avoid an excessively long simulation runtime. In other words, we assume that the system has a perfect correction for artifacts due to crosstalk and pileup.

Simulating voxelized phantoms in CatSim requires files that describe the phantom in terms of resolution and size while referring to material density matrices that provide the fractions of a given material at each image point. A MATLAB script was used to convert the raw XCAT output to the format required by CatSim. To reduce memory and compute time requirements, a material decomposition into polyethylene (PE) and polyvinyl chloride (PVC) was performed to reduce the number of density matrices to two. The choice of basis materials was motivated by PE and PVC being commonly used calibration materials since they can be linearly combined to approximate human tissues with positive coefficients.³²

The linear attenuation coefficient at any point in the body can be approximated by a linear combination of basis functions, which can be expressed as

The material decomposition into PE and PVC was performed by first calculating their respective linear attenuation coefficients at 40 and 70 keV using the “GetMu” function available in the XCIST repository on GitHub, as well as corresponding coefficients for the materials in the XCAT phantom. For each point, an equation system

The forward projection was simulated in CatSim, and the image reconstruction was performed using a proprietary reconstruction prototype software for photon-counting CT. Based on registered photons in each energy bin, a material decomposition into PE and PVC was computed using a projection-space maximum-likelihood-based material decomposition algorithm, and basis images were reconstructed through filtered back-projection with a standard kernel. The outputs of the reconstruction software were water and iodine basis images, i.e., by default the software performed a change of basis from PE and PVC to water and iodine. This step was reversed using the inverse basis transformation matrix, and we used the initial PE and PVC basis images to compute the final reconstructed VMIs. The reconstruction field of view was set to 350 mm and the reconstructed slice thickness was 0.4167 mm. We did not apply any noise reduction algorithm to the image data.

Inspired by the work of Näsmark and Andersson, as described in Ref. 33, VMIs of 40 and 70 keV were used as input to the network. The basis images were therefore converted to 40 and 70 keV VMIs through linear combination, consistent with Eq. (6).

2.4.

Neural Network Model

The network used in this project was a U-Net inspired by a network architecture employed as a component of denoising diffusion probabilistic models (DDPMs), as described by Song et al. in Ref. 34, who based their architecture on Ho et al. in Ref. 35. While using this network as a starting point, adaptations were made to employ the architecture for image-to-image translation that could be trained with the $L1$, MSE, and VGG16 loss functions. A necessary change was to remove the time embeddings used for the score-matching in diffusion models. The network had one BigGAN residual block per resolution and did not employ the anti-aliasing based on finite impulse response for up- and down-sampling, as was proposed in Ref. 34.

The ADAM optimizer was used with a learning rate of ${10}^{-4}$, ${\beta }_1=0.5$, and ${\beta }_2=0.9$. Common loss functions used for neural network applications to CT image data, such as for denoising, are explored in Ref. 36. To evaluate which loss function was most suitable for this task, four loss functions were evaluated after 100 epochs: MSE loss, $L1$ loss, VGG16 loss with feature extraction at the ninth layer, and a weighted combination of VGG16 and $L1$ loss, hereafter referred to as VGG16_L1 loss. Because VGG16 is trained with RGB images, the SPR map was repeated to three channels and normalized to match the statistics of the VGG16 training data, before being given as input to VGG16. The final network was trained with the VGG16_L1 loss, defined as

We used two input channels because VMIs at two distinct energies were provided as input, and a single output channel because the corresponding ground truth SPR map is a single grayscale image. The original images were of $1024\times 1024$ resolution, which required more GPU memory than what was available. Images were therefore re-sized to a resolution of $256\times 256$. The network was subsequently trained on a Razer Blade 14, with an AMD Ryzen 9 5900 HX (3.30 GHz) processor with an NVIDIA GeForce RTX 3080 graphics card, and the training took $\sim 12\mathrm{h}$.

Because of the stochastic gradient descent with random draws from the examples, as well as random initializations of the network parameters, every training instance results in a slightly different model. To get a sense of this uncertainty, the training was repeated fifteen times with identical data sets, settings, and noise realizations.

The XCAT phantom was transformed into a synthetic 70 keV VMI, using the theoretical attenuation coefficients computed from each material’s elemental composition. The attenuation map was then transformed into a synthetic VMI at 70 keV using the relationship

During the simulation process, the phantom is rescaled in relation to the original. To ensure pixel-to-pixel correspondence between the simulated CT image and SPR map prior to training, image registration was performed in MATLAB. In MATLAB, the built-in functions “imregconfig” and “imregtform” were used to generate a registration transform, fitting the synthetic VMI to the corresponding simulated VMI. The resulting registration transform was then applied directly to the SPR maps because these were of identical scale as the synthetic VMI generated from the XCAT phantom.

The data set was randomized and divided into a training and test set with a ratio of 80 to 20, using the built-in function “train_test_split()” from scikit-learn in Python. The training set was then divided into a training and validation set with the same ratio.

2.5.

Non-Deep Learning Methods

To compare the accuracy of the network to that of conventional methods, a stoichiometric calibration and a VMI-based method originally adapted for DECT (hereafter referred to as N&A) were performed on the simulated PCCT data. A detailed description of these implementations is provided in the Appendix.

3.1.

Simulations

To evaluate the accuracy of the photon-counting simulation, a synthetic VMI at 70 keV was computed from an XCAT phantom based on the material in each voxel. The HUs of a row through the image were plotted with the HUs of the corresponding row and slice from the simulated 70 keV VMIs. The results are shown in Fig. 2.

Fig. 2

Plotted HUs for a true 70 keV VMI based on the XCAT phantom and a simulated 70 keV VMI from CatSim (WL: 0, WW: 2000).

The RMSE of $\frac{\mu }{{\mu }_{\text{water}}}$ was calculated for each ROI, with $\mu$ calculated from the HU in the simulated VMI by inverting Eq. (7). The results are shown in Table 2.

Table 2

RMSE of μμwater calculated for each material’s ROIs.

3.2.

Performance of the Neural Network

The four different loss functions were each used to train the network and were subsequently evaluated. The mean RMSE for three regions of interest was calculated for the whole brain material for each loss function. The results are shown in Table 3.

Table 3

RMSE for the whole brain material for each of the loss functions.

For further evaluation, the results of the best-performing loss function VGG16_L1 were used. Results for observed VMI, predicted SPR, true SPR, and their difference are shown for one selected slice in Fig. 3.

Fig. 3

Observed 40 keV VMI (a) and 70 keV VMI (b), predicted (c), and true (d) SPR maps and a difference image (e) where true SPR values have been subtracted from the predicted values.

To evaluate the SPR prediction, a slice and a row were selected and the predicted and true SPR values were plotted. The result is shown in Fig. 4.

Fig. 4

Comparison of predicted and true SPR values.

Errors in SPR were calculated for tissues of lower density (HU below 0), higher density (HU above 100), and tissues in the range of 0 to 100 HU. Because most tissues fell into the 0 to 100 HU range, these were grouped into intervals of 10 HU for better visualization. Histograms of the error distributions are shown in Fig. 5.

Fig. 5

SPR error distributions for different HU intervals from the 40 keV (a) and 70 keV (b) VMIs. The parentheses denote the number of data points in the given interval. The maximal error was 0.59 and the minimal error −0.68.

The slice from Figs. 3 and 4 was evaluated over fifteen separate training sessions. Errors in SPR were calculated for white matter, gray matter, the whole brain material, and the skull, in three ROIs each. The RMSEs are shown in Table 4. For the relative error calculation, the mean predicted SPR in the ROI was taken as the predicted value. The results are shown in Table 5. The relative standard deviations, referring to a sample standard deviation for the ROI divided by the true SPR value, were calculated for each ROI and are shown in Table 6.

Table 4

Mean RMSE and standard deviation for ROIs of the white matter, gray matter, whole brain, and skull materials, in the SPR map predicted by the neural network. The column “mean” refers to the mean RMSE calculated over the three ROIs for one material.

Table 5

Mean relative errors and standard deviations for ROIs of the white matter, gray matter, whole brain, and skull materials, in the SPR map predicted by the neural network. The column “mean” refers to the mean relative error calculated over the three ROIs for one material.

Table 6

Mean and sample standard deviation of the relative standard deviations for ROIs of the white matter, gray matter, whole brain, and skull materials, in the SPR map predicted by the neural network. The column “mean” refers to the mean relative standard deviation calculated over the three ROIs for one material.

Errors for the SPR maps estimated from the stoichiometric calibration and the N&A method were evaluated for the same slice and ROIs. The RMSEs are shown in Table 7, relative errors in Table 8, and relative standard deviations in Table 9. The optimal VMI pairs for the N&A method were 46 and 55 keV for the soft tissue and 51 and 61 keV for the bone.

Table 7

RMSE for ROIs using SPR calculation with the N&A method and SECT stoichiometric calibration.

Table 8

Relative errors for ROIs using SPR calculation with the N&A method and SECT stoichiometric calibration.

Table 9

Relative standard deviations for ROIs using SPR calculation with the N&A method and SECT stoichiometric calibration.

6.1.

Stoichiometric Calibration

It has been reported that a lack of consensus regarding the implementation of stoichiometric calibration has led to treatment center inter-variability.³ To avoid ambiguity, the standardized guide by Peters et al. in Ref. 37 was closely followed, using the available code and Excel input template on GitHub ( https://github.com/CTinRT/HLUT-guide). A Gammex Advanced Electron Density phantom was defined in CatSim, with insert material compositions in accordance with those listed in the input template on GitHub. The elemental ionization energies were adjusted manually to correspond to the values used for the other SPR calculations in this work.

PCCT acquisitions to specify the HU for lung, soft tissue, and bone inserts were simulated in CatSim for several insert configurations, as suggested in Ref. 37. The CT protocol parameters were the same as for the XCAT CT acquisition. Twelve slices were acquired for each configuration. Reconstructed basis images for PE and PVC were then recombined into 70 keV VMIs, and the HU for a tissue insert was calculated as the mean value of a circular ROI covering ca. 70% of the insert’s inner cross-sectional area, averaged over the 12 slices.

The stoichiometric calibration was performed with the code available on GitHub. The only changes made to the script were the Mayneord factor, used for the effective atomic number calculation, which was set to 3.21 (as in Ref. 33), and $\beta$ which was set manually to the same value as for the other methods in this work. The HLUTs were then calculated separately for the head and body phantom. Because this work only included XCAT phantoms of the head, the head HLUT was selected for further evaluation. Based on the given breakpoints of the line segments, a piecewise linear function was defined using the “pwlf” package in Python. The resulting function was subsequently used to predict SPR maps based on simulated 70 keV VMIs of the XCAT phantom.

6.2.

SPR Estimation from Pairs of VMIs

A VMI-based method originally adapted for DECT was applied to the simulated PCCT image data. The proposed method by Näsmark and Andersson, described in full detail in Ref. 33, employs pairs of VMIs to calculate the effective atomic number and relative electron density. The mean ionization energy is derived from the effective atomic number via separate parameterizations for soft and bone tissues, after which the mean ionization energy and relative electron density are used as input quantities to Eq. (3).

In the same work, the authors demonstrated that the SPR accuracy was highly dependent on the VMI pair, and the optimal choice depends on the scanner’s ability to generate accurate VMIs for different energies. An initial optimization was performed on VMI pairs in the range of 40 to 140 keV, where each pair was assessed based on its accuracy in estimating the SPR of the tissue inserts. The optimal pair for each tissue category (lung tissue, soft tissue, and bone) was selected as the one yielding the smallest RMSE, as compared with the ground truth calculated with Eq. (3) based on the elemental composition of the inserts and listed reference tissues from Ref. 26. We did not have to explicitly enforce the skewness criterion recommended in Ref. 39 because this was fulfilled by the highest ranking VMI pair anyway. The same Gammex phantom image data obtained for the stoichiometric calibration was used for the optimization. Because this work included only XCAT head phantoms, the optimization was based on the Gammex head phantom acquisitions.

6.3.

Effects of Crosstalk and Pileup

We chose to exclude crosstalk and pileup in the simulations, which are known to deteriorate the image quality in PCCT. This choice was motivated by the fact that the effects are negligible in the final reconstructed image, after the application of correction algorithms in the image chain. We simulated a scan of a water phantom to validate this assumption, by first excluding crosstalk and pileup, and then including these effects. A pileup correction algorithm was applied in the reconstruction.⁴⁰ The example is shown in Fig. 6.

Fig. 6

Reconstructed VMIs at 70 keV for a simulation excluding (a) and including (b) crosstalk and pileup, as well as line profiles through the widest part of the cylinder for the simulation excluding (c) and including (d) crosstalk and pileup.