Development of a Monte Carlo-based scatter correction method for total-body PET using the uEXPLORER PET/CT scanner

The development and evaluation of the SC algorithm were performed using the uEXPLORER total-body PET scanner, which comprises eight rings (also called units), each with an AFOV of approximately 24 cm. Each unit incorporates 24 transaxial detector modules, comprising of 14 (axial) by 5 (transaxial) detector blocks each of which consists of 6 (axial) by 7 (transaxial) scintillator crystals. The crystals are made of LYSO and have dimensions of 2.76 × 2.76 × 18.1 mm³, with a pitch of 2.85 mm (Omidvari et al 2022). The scanner records coincidences with a maximum unit difference (MUD) of 4, resulting in a maximum acceptance angle of 57 degrees. The coincidence time window is variable, ranging from 4.5 (MUD = 0, intra-unit coincidences) to 6.9 ns (MUD = 4, full acceptance angle). The detectors have an energy resolution of 11.7% and the system uses an energy window of 430–645 keV for coincidence detection. Data are collected in list-mode format with a TOF resolution of 505 ps.

The SC algorithm was developed and implemented in an iterative PET image reconstruction framework. This in-house tool was generated in collaboration with the scanner manufacturer to allow direct incorporation of raw list-mode data acquired from the uEXPLORER scanner. Our tool allows the use of independently generated data corrections as well as full adaptation of reconstruction parameters outside of the clinical environment. The total-body PET reconstruction environment utilizes TOF list-mode ordered-subsets expectation maximization (LM-OSEM) (Leung et al 2021a) the implementation of which is based on Reader (1998). More details on the implementation are given in the supplemental materials.

The goal of this study was the development of an in-house MC-based method to estimate the number of scattered events contributing to each LOR and TOF bin. The framework is outlined in the flow chart in figure 1; the SC estimation was performed between individual reconstruction iterations. Thus, in the first step of the reconstruction framework, an initial image estimate was created using one OSEM iteration with 13 subsets, considering corrections for random events, attenuation, and dead time, but excluding scatter. Only the raw list-mode data and vendor-provided normalization factors were used as input, and all listed corrections have been developed by the UC Davis research team. This estimate served as the activity distribution map for the MC simulation toolkit SimSET (version 2.9.2—customized for this application). The spatial distribution of true and scattered events measured by the PET scanner was then simulated, based on the initial image estimate and an attenuation map. The latter was derived from the Hounsfield units of a CT image that were transformed into attenuation coefficients using bilinear interpolation (Carney et al 2006). The recorded events were binned into block-based sinograms, with each block comprising 6 axial and 7 transaxial crystals.

Zoom In Zoom Out Reset image size

To correctly estimate the amount of scatter contributing to each detector block pair, the simulated events were scaled to the measured events. Subsequently, scatter estimates for all list-mode events were calculated and used in the next OSEM iteration update. That way the subsequent image reconstruction iteration will result in a first scatter corrected image estimate, which in turn will serve as next updated image guess for the scatter simulation. Since the initial image estimate—that is used as input for the simulation—contains scatter, an initial over-correction of scatter in the first SC iteration is expected as described by Levin et al (1995). Consequently, it is plausible that the second SC iteration may lead to an underestimation of the scatter contribution in the PET image. Therefore, the SC process was repeated three times, refining the correction terms in each step. After the completion of SC factor calculation, the OSEM reconstruction was performed with all corrections included (details below).

The MC simulation toolkit SimSET requires a set of input parameters, including a three-dimensional (3D) attenuation map and an activity distribution map, to simulate the detector response. Both maps had an image size of 239 × 239 × 679 voxels, with an isotropic voxel size of 2.85 × 2.85 × 2.85 mm³. To preserve the granularity of the continuous attenuation map generated by CT, all materials were defined through densities ranging from 0.01 to 2.2 g cm⁻³ (in increments of 0.01 g cm⁻³) instead of using the 29 standard materials defined in SimSET. Water was chosen as material for all voxels and was assigned a density that had the closest corresponding 511 keV attenuation coefficient. Then, the water density-map was converted into the material index map required by SimSET (Leung 2021).

In the SimSET toolkit the geometry of the uEXPLORER scanner was defined as described above. Photons with an energy of 511 keV were created for each decay, and the detector energy resolution was defined as 11.7% (Zhang et al 2020). In the simulation, intra- and inter-crystal scatter was simulated, however, no dead time effects were considered. The uEXPLORER does not have septa between blocks or crystals, nor end shields that are used in shorter scanners to limit scatter contributions from outside the FOV. Therefore, no hardware components had to be simulated.

In each iteration, 2.5 billion decays of ¹⁸F were simulated—distributed over 44 instances executed on two Intel^® Xeon^® Gold 6126 CPUs at 2.6 GHz in parallel. Each SimSET instance collected data from about 56.8 million decays. No variance reduction methods such as stratification or forced detection were used and the total computation time for the simulation was ~1 h for each iteration, making it a total of 3 h of computation time for SC. The SimSET output files with scattered and true events were converted into simple list-mode files, where each recorded event was a 10-byte binary string encoding the axial and transaxial coordinates of the two coincident crystals, as well as the TOF information (2 bytes for each information). During that process, a random number-based Gaussian TOF blurring was applied to each event to simulate the uEXPLORER's TOF resolution of 505 ps (Spencer et al 2020).

The sinograms of coincident events obtained from a high-resolution total-body PET scanner at the crystal level are too large to manage due to large number of crystals, large acceptance angle, and enormous number of approximately 90 billion possible LORs. At typical scan conditions, crystal-based sinograms would provide low statistical significance and would exhibit substantial noise. Therefore, a strategy was devised to overcome this in which the recorded events were binned into block sinograms corresponding to the detector blocks in the actual scanner consisting of 6 × 7 crystals. A five-dimensional sinogram was created with 91 radial and 60 angular bins per detector plane, 112 × 112 planes in total, and 27 TOF bins. The sinogram TOF bins combined seven TOF bins from the measured raw data with 39.0625 ps bin width, therefore having a width of 273.7 ps, which corresponds to approximately half the scanner's coincidence time resolution. Despite the use of small 16-bit encoding for the bin entries, the file size of each 5D sinogram was approximately 3.5 GB, due to the enormous number of projections employed in its generation.

The measured PET projection data, or prompt events (P), consist of true (T), scattered (S), and random (D) coincidences, the latter of which are estimated using delayed coincidence channels. Based on the activity distribution in the FOV, the scattered and true events are simulated, the distribution of which is expected to approximate that of the measured prompt minus random events. Since the simulation only provides count distributions relative to the total number of simulated events, correct scaling of the scatter sinograms is required to obtain the absolute contribution of scattered events to each sinogram block in the measured data. To that end, a scaling factor ${f}_{k}$ for each sinogram k (with k = 1, 2, ... 112 × 112) must be introduced. Scatter scaling follows the assumption of approximately equal SF in simulated and measured sinograms:

$$\begin{eqnarray}&&\begin{array}{c}{\mathrm{SF}}_{k}=\frac{\displaystyle {\sum }_{b}{S}_{b}^{\mathrm{meas}}}{\displaystyle {\sum }_{b}{\left({P}_{b}-{D}_{b}\right)}^{\mathrm{meas}}}=\frac{\displaystyle {\sum }_{b}{S}_{b}^{\mathrm{sim}}}{\displaystyle {\sum }_{b}{\left({T}_{b}-{S}_{b}\right)}^{\mathrm{sim}}}\\ \iff \displaystyle \sum _{b}{S}_{b}^{\mathrm{meas}}=\displaystyle \sum _{b}{S}_{b}^{\mathrm{meas}}\cdot \frac{\displaystyle {\sum }_{b}{\left({P}_{b}-{D}_{b}\right)}^{\mathrm{meas}}}{\displaystyle {\sum }_{b}{\left({T}_{b}-{S}_{b}\right)}^{\mathrm{sim}}}=\displaystyle \sum _{b}{S}_{b}^{\mathrm{meas}}\cdot {f}_{k},\end{array}\end{eqnarray} \tag{ 1 }$$

where b denotes the block sinogram bin number (with b = 1, 2, ... 91 × 60). Since the distribution of scattered and true events over the individual TOF bins is not expected to be different for simulated and measured data, the entries from all TOF bins have been summed together for each block bin b in equation (1) to reduce noise, and a single scaling factor for all TOF bins has been used. The scaling factor for each plane k is therefore the number of measured prompts minus delayed events divided by the number of simulated true and scattered events that are contained in their corresponding sinograms. Here, all negative entries in the P-minus-D sinogram were set to zero prior to scaling. No other sinogram smoothing was applied besides the block-based binning. The same scaling factor is used for all TOF bins of one sinogram plane, and the number of scattered events in the measurement contributing to sinogram bin b and TOF bin t can then be estimated using scaled simulated scattered events:

$$\begin{eqnarray}&&{S}_{b,t}^{\mathrm{meas}}\approx {S}_{b,t}^{\mathrm{sim}}\bullet {f}_{k}.\end{eqnarray} \tag{ 2 }$$

Figure 1 in the supplementary materials plots the average scaling factor versus axial block ring difference for the phantom and the human subject study discussed in this work. The strong fluctuations emphasize the need for calculating a scaling factor for each individual block plane instead of using a global scaling factor for the whole data set.

Based on the size and activity distribution within the field of view (FOV), along with the object's material composition, it is noteworthy that only a fraction (10%–30%) of all simulated decays results in the recording of true or scattered coincident events. To limit the amount of noise introduced by the scatter estimate into the next OSEM iteration, Gaussian smoothing was applied to the sinograms of scattered events using a smoothing kernel with a full width at half-maximum (FWHM) of 8 bins. Since this step was performed after the sinogram scaling, the preservation of the total number of events in each sinogram plane was ensured.

The average number of scattered events contributing to LOR i located in detector block pair b with TOF bin t is approximated following the assumption that there is a uniform distribution of scatter over all LORs in one block sinogram bin:

$$\begin{eqnarray}&&{s}_{i,t}\approx \frac{{S}_{b,t}^{\mathrm{sim},\mathrm{Gauss}}\bullet {f}_{k}}{{N}_{\mathrm{LOR}}},\end{eqnarray} \tag{ 3 }$$

where ${S}_{b,t}^{\mathrm{sim},\mathrm{Gauss}}$ is the number of entries in TOF bin t and block sinogram bin b in the smoothed 5D sinogram of simulated scattered events, and ${N}_{\mathrm{LOR}}={\left(6\times 7\right)}^{2}=\mathrm{1,764}$ denotes the number of LORs between two detector blocks. The obtained values from equation (3) served as SC terms in the TOF LM-OSEM algorithm.

Table 1 in the supplementary materials provides an overview of all the parameters used in the SC algorithm described above.

First, three NEMA image quality (IQ) phantoms and a uniform cylinder phantom with a 20 cm diameter and a length of 32 cm were placed on the scanner bed as shown in figure 2(A). The configuration is shown in the CT scan in figure 2(A); the phantoms are scanned on the uEXPLORER total-body PET/CT scanner simultaneously in one bed position. The phantoms are labeled with numbers 1 (top end) to 3 (axial center) to indicate their axial position in the later analysis. Following the NEMA NU 2–2018 protocol (NEMA 2018), the background concentration in all three NEMA IQ phantoms and in the homogenous cylinder was 5.24 kBq ml^–1. The activity concentration ratio between the hot spheres of the IQ phantom and the background was 3.9:1. The lung insert in the first IQ phantom was filled with Styrofoam balls and air whereas the other two phantoms were filled with a mixture of Styrofoam balls and water to investigate the scatter removal in cold regions for various attenuation levels.

Zoom In Zoom Out Reset image size

Additionally, a 79 year old male with non-specific chronic granulomatous disease and incidentally discovered pulmonary nodules was injected with 377 MBq (10.2 mCi) [¹⁸F]FDG and scanned for 20 min at 2 h post injection following the currently used clinical scan protocol (Nardo et al 2020, 2021). To provide a more challenging test for the SC algorithm, the chosen subject had an above-average body mass index (BMI) of 32.5 kg m⁻², where higher scatter fraction and prominent photon-attenuation can lead to a degradation of image quality. The participant's data were acquired after informed written consent following IRB approval (IRB #1506448).

Images were reconstructed using 4 OSEM iterations with 13 subsets with one iteration of MC-based scatter estimate in between OSEM iterations to update the SC terms. The data was reconstructed into an image matrix with a size of 239 × 239 × 679 voxels with an isotropic voxel size of 2.85 mm. No point spread function (PSF) modeling was applied. For comparison, all reconstructions were also performed without SC. Images of the human subject were smoothed with a Gaussian filter with a FWHM of 4.0 mm.

All phantom images were visually inspected for scatter removal and artifacts related to SC and scatter corrected images were compared to their non-scatter-corrected counterparts. The scatter removal in the phantom's plastic walls and inside the cold lung insert were investigated first visually and then with a line profile in the transverse plane where the spheres were located.

The quantitative bias of the activity concentration ratio between the largest hot sphere and the warm background was calculated. A spherical region-of-interest (ROI) with a diameter of 30 mm was drawn in the 37 mm sphere of each IQ phantom and three axially oriented cylindrical ROIs with a diameter of 50 mm and a length of 150 mm were drawn in the warm background (BG). This quantitative analysis method is different from the calculation of contrast recovery coefficient (CRC) (e.g. following NEMA NU 2–2018 (NEMA 2018)). in the sense that it solely estimates the bias on the activity concentration ratio and therefore serves as a measure of quantitative accuracy of corrections. Furthermore, since an ROI diameter smaller than the sphere was chosen, the calculation does not suffer from potential partial volume effects at the boundary of the sphere (i.e. dependence on voxel size). The ratio of the mean number of counts in the hot ROI and in the two ROIs in the BG was calculated and the relative deviation to the actual activity concentration ratio in the phantom was calculated.

To investigate the axial uniformity of the corrections, the bias on the background activity concentration ${c}_{\mathrm{IQ}}\,$ of each individual NEMA IQ phantom was calculated with respect to the average activity concentration in the uniform cylinder phantom ${c}_{\mathrm{cyl}}:$

$$\begin{eqnarray}&&\mathrm{bias}({\mathrm{BG}}_{\mathrm{IQ}})=\frac{{c}_{\mathrm{IQ}}-{c}_{\mathrm{cyl}}}{{c}_{\mathrm{cyl}}}.\end{eqnarray} \tag{ 4 }$$

A cylindrical ROI with a diameter of 120 mm and a length of 220 mm was drawn in the homogenous cylinder phantom and the mean value was calculated. The accuracy of the SC should be independent of axial position and therefore the bias should be consistently minimal for all three phantoms.

Following NEMA NU 2–2018 (NEMA 2018) the CRC was calculated for all three phantoms as well as the residual error in the lung insert. Both evaluations were performed for scatter corrected and non-scatter-corrected reconstructed images.

Human subject images were visually inspected for scatter removal and artifacts related to SC. In particular, the region around the bladder was examined where over-correction of scatter frequently occurs in clinical practice, potentially obliterating tumor uptake in pelvic-located tumors as described by Roman-Jimenez et al (2016). Maximum-intensity-projections (MIPs) were created for both scatter corrected and non-scatter-corrected reconstructed images. Additionally, line profiles through image slices were examined at three distinct locations: (i) in the sagittal view we investigated scatter removal in the trachea; in the coronal view we investigated scatter removal (ii) in the air between the neck and arms, (iii) between the thighs, and (iv) inside a metal prosthesis in the patient's right knee. We quantified the scatter clearance by identifying the minimum activity concentration in the four line profiles on both scatter corrected and non-scatter-corrected images and calculated the relative difference between the two with respect to the non-scatter-corrected image.

Figure 2 shows the CT image (A) and the UCD in-house reconstructed PET image (B) of the three NEMA IQ phantoms and the homogenous cylinder. Figure 2(C) shows a zoomed view of the second IQ phantom. No residual scatter in the plastic parts of the phantom is visible. The two red arrowheads point to the residual activity in the spheres' filled stems, which has correctly been reconstructed showing activity, while any scatter in the surrounding plastic has been removed successfully. The dashed line in subfigure(C) indicates the line profile that is shown in subfigure(D). Since scatter is broadly distributed in the transaxial direction and can in some cases be approximated with a Gaussian distribution, the non-scatter-corrected image shows a scatter-induced gradient towards the center of the transverse FOV. In contrast, the line profile on the scatter corrected image is expectedly flatter. The minimum activity concentration in the line profile inside the lung insert was 1.102 kBq ml⁻¹ at a transaxial position of 352 mm for the non-scatter-corrected image, and 0.059 kBq ml⁻¹ at 346 mm for the scatter corrected image, respectively. This constitutes a reduction of the minimum concentration by 94.6% and brings the voxel values inside the lung insert close to the ground truth of 0.000 kBq ml⁻¹.

Figure 3 shows the bias on the activity concentration ratio between the largest hot sphere and the warm background ROIs and the bias on the background concentration with respect to the homogenous cylinder phantom. The bias on the concentration ratio was between −1.35% for the NEMA IQ phantom at the head end of the FOV and −3.15% for the phantom in the axial center. The bias on the background concentration varied between −0.75% for the phantom at position two and +0.63% for the phantom in the axial center.

Zoom In Zoom Out Reset image size

Figure 4 compares the CRC and background variability measured for all three NEMA IQ phantoms, with and without SC. For all three phantoms the CRC for all sphere sizes is higher on the scatter corrected images compared to the non-scatter-corrected ones. For the 10 mm sphere the measured CRC for the three NEMA phantoms was between 52.9% and 55.7% with SC and between 46.7% and 50.0% without SC, respectively. For the 37 mm sphere the CRC values were between 83.7% and 86.0% with SC and between 74.1% and 77.1% without SC. The average relative CRC improvements with SC compared to the non-scatter-corrected images were 17.2% for phantom (1), 7.2% for phantom (2) and 11.1% for phantom (3), respectively.

Zoom In Zoom Out Reset image size

The background variability in all phantoms is consistently improved when SC is applied—except for one data point observed in phantom (1). Here, on the scatter corrected images an average relative reduction of the background variability by 10.2%, 18.7% and 34.3% were obtained for phantoms (1), (2) and (3), respectively.

Figure 5 shows the residual error in the lung insert in the NEMA IQ phantoms at all three axial positions. While without SC the residual error was between 10.43% for the phantom at position one and 16.91% for position two, the values with SC ranged between 1.26% for phantom 1 and 2.08% for phantom 3. This implies a relative reduction of the residual error between 87.9% and 91.4% when SC was applied with respect to the non-scatter-corrected case.

Zoom In Zoom Out Reset image size

Figure 6 shows reconstructed images of the human subject with (A) and without (B) scatter correction. Qualitatively, the scatter corrected image shows improved contrast and artifact-free correction of scatter around the urinary bladder. Subfigures (C) through (F) show the position of line profiles drawn on single image slices through the trachea in sagittal view (C), through the neck and arm region (D), through the thighs (E) and through the prosthesis in the knee (F) in coronal view.

Zoom In Zoom Out Reset image size

In the first line profile the trachea was located at a sagittal position between 382 and 392 mm. The minimum activity concentration in the line profile through the trachea was 0.046 kBq ml⁻¹ in the scatter corrected image and 0.350 kBq ml⁻¹ in the non-scatter-corrected image, respectively, which constitutes a reduction of the minimum activity concentration in the trachea by 86.9%. In the line profile through the neck and arm region, the residual activity concentration was evaluated at two different locations: the transaxial position of the air gap between head and right arm was located between 284 mm and 303 mm. The line profile showed a minimum activity concentration of 0.018 kBq ml⁻¹ with SC and 0.259 kBq ml⁻¹ without SC, respectively, which is a reduction by 93.1%. The gap between head and left arm was found between 435 mm and 463 mm and the minimum activity concentration was 0.002 kBq ml⁻¹ with SC and 0.153 kBq ml⁻¹ without SC, respectively, which implies a reduction by 98.7%. Finally, the air gap between the thighs was found at a transaxial position between 360 mm and 380 mm and a minimum activity concentration of 0.003 kBq ml⁻¹ with SC and 0.150 kBq ml⁻¹ without SC was obtained, which constitutes a reduction by 98.0%.

The amount of scatter that occurred inside and around the right knee prosthesis was noticeably reduced (see red arrowheads in figure 6) and residual uptake around the prosthesis was found to have the expected level for such an intervention. The line profile through the prosthesis (F) showed reduced residual activity in the femoral component of the prosthesis indicated by the black arrow in the lateral and by the red arrowhead in the medial component, respectively. Especially the medial aspect exhibited a reduction of activity concentration from a peak value of 8.38 kBq ml⁻¹ in the non-scatter-corrected image down to 0.57 kBq ml⁻¹ in the corrected image at the same transaxial position. This constitutes a reduction of 93.2%.

Importantly, no washout effects were found surrounding areas with high uptake, and in all the line profile evaluations the minimum activity concentration was greater than 0.000 kBq ml^–1 after SC. This indicates that the SC algorithm did not overcorrect for scatter in these regions. A visualization of the scatter contribution to the different body regions in the reconstructed image can be found in the supplemental materials, figure 2.