The Mediso AnyScan Trio SCP (SPECT-CT-PET) system is a hybrid imaging system with three imaging modalities. The SPECT component of the model installed at NPL has a NaI(Tl) scintillation crystal thickness of 9.5 mm and three detector heads. It can operate either in triple- or dual-head acquisition mode where the detector heads are at a 120\(^\circ\) or 180\(^\circ\) separation, respectively. In dual-head acquisition mode, the third head is placed at 90\(^\circ\) to the other two, but does not acquire counts. A photograph of the system in dual-head acquisition mode is shown in Fig. 1. Two collimators were considered in this work: Low Energy High Resolution (LEHR) and Medium-Low Energy General Purpose (MLEGP) collimators, recommended for imaging with \(^}\)Tc and \(^\)Lu, respectively. A description of their hole sizes and septal thicknesses from the Mediso AnyScan catalogue [15] is given in Table 1.

Table 1 The specifications of the low energy high-resolution (LEHR) and medium-low energy general purpose (MLEGP) collimators used in this workExperimental data acquisitionPhantom acquisitions

Experimental acquisitions of \(^}\)Tc and \(^\)Lu were performed with LEHR and MLEGP collimators, respectively. Images of each isotope were acquired with a uniform and a localised distribution of activity. A cylinder with 20 cm internal diameter and height was filled with a uniform solution of \(^}\)Tc. A commercial cylindrical Jaszczak phantom was used for the \(^\)Lu acquisition; with an internal height of 18.6 cm and a diameter of 21.6 cm [16]. A commercial NEMA phantom [17] with six active spheres in an inactive background of water was used for both isotopes. The same activity concentration was dispensed to all spheres. Water was used as a scattering medium due to its similar density to human tissue. All phantom scans were performed in dual-head acquisition mode with a fixed detector radius of 250 mm; they are summarised in Table 2.

One emission (EM) and one scatter (SC) windows were simultaneously acquired for \(^}\)Tc imaging. A larger (5%) scatter window was used for the cylindrical phantom. 20% emission and 6% scatter windows were acquired for the two emissions of \(^\)Lu [18, 19]. The Mediso AnyScan system permits a maximum of four windows per scan so two consecutive scans with one emission and the adjacent scatter windows were performed.

Table 2 Details of the phantom scans performed

Source activities were measured with the NPL secondary standard ionisation chamber system (Vinten 671). The uncertainty budgets for the NEMA phantom activities are given in Table 3. A final uncertainty of 0.892% and 0.503% was obtained, through adding the components in quadrature, for \(^}\)Tc and \(^\)Lu, respectively.

Table 3 The percentage uncertainty components for activity per unit mass of \(^}\)Tc and \(^\)Lu solutions used in the NEMA phantom

All experimental images were acquired in DICOM format with 128\(\times\)128 matrices and voxel sizes of (4.2578 mm)\(^3\). A CT scan was performed directly following each SPECT scan, with 512\(\times\)512 matrices of \(0.9766 \times 0.9766 \times 0.625\) mm\(^3\) voxels.

Extrinsic spatial resolution

Point source measurements of \(^}\)Tc were obtained to measure the extrinsic spatial resolution of the system with LEHR collimators. An active volume of approximately 0.5 ml and 12 MBq was drawn into a 2.5 mL syringe which was positioned in a central, reproducible position on a custom-made perspex board. A single detector head was positioned directly above the patient bed, at the minimum height permitted by the system (106.5 mm). A static planar acquisition was conducted with a termination condition of 500,000 counts in an energy window of 140.5 keV ± 10%. This process was repeated for three additional head radii: 150, 250 and 350 mm. Further measurements were made with the head at 106.5 mm height and the source shifted laterally (150 mm in the long axis of SPECT head - ‘CentreLeft’ and −150 mm in the long axis and 100 mm in the short axis - ‘TopRight’). All experimental acquisitions used a 1024\(\times\)1024 matrix with pixel width 0.5322 mm to reduce the uncertainty on experimental position in the scan (half the acquisition pixel size, or 0.2661 mm). The experimental acquisitions were then down-sampled to 128\(\times\)128 matrices to match SPECT pixel size used in the rest of this work.

Background radiation

Extrinsic background acquisitions with LEHR and MLEGP collimators were performed to record counts with no source present. Tomographic acquisitions were acquired for all energy windows used in the phantom acquisitions for 32 s per projection. The energy spectra were extracted directly using manufacturer-provided software and are shown in Fig. 2.

Intrinsic camera efficiency

Standard calibration sources were used to determine the experimental intrinsic efficiency of the Mediso SPECT AnyScan system, these were \(^\)Cd, \(^\)Ba, \(^\)Ce, \(^\)Pb, \(^\)Am, \(^}\)Tc, \(^}\)Te, and \(^\)Lu. A liquid solution of each radionuclide was dispensed to a 2 ml ISO ampoule [20] and placed in a reproducible position with the SPECT head at a fixed radius of 455 mm. No collimators were used. All calibration source data were acquired with a \(^}\)Tc energy-correction map [21]. A photograph of the set-up is shown in Fig. 3. All sources were traceable to national primary activity standards, with activities as shown in Table 4. Static acquisitions were acquired for 600 s for each source, apart from \(^\)Cd and \(^\)Pb which were acquired for 3600 and 7200 s, respectively, to account for their low count rates.

Table 4 The \(\upgamma\) emissions used to calculate the intrinsic detector efficiency

A Gaussian fit was applied to each photopeak in the energy spectrum of each isotope. An absolute detection efficiency was then calculated for each of the peaks according to

$$\begin \epsilon = \frac \end$$

(1)

where N is the net count rate in the peak, A is the source activity, and \(I_\gamma\) is the emission intensity of the \(\upgamma\) ray of interest. The uncertainty on the absolute efficiency, \(\sigma _\epsilon\), was determined through standard uncertainty propagation, assuming uncorrelated variables, as

$$\begin \sigma _\epsilon ^2 = \left| \frac\right| ^2 \sigma _N^2 + \left| \frac\right| ^2 \sigma _A^2 + \left| \frac }\right| ^2 \sigma _}^2 = \epsilon ^2 \left[ \left( \frac}\right) ^2 + \left( ~\frac}\right) ^2 + \left( \frac}}}\right) ^2 ~\right] \end$$

(2)

where \(\sigma _N\), \(\sigma _A\) and \(\sigma _}\) are the uncertainties on the net count rate, activity and emission intensity, respectively.

Simulation model

A full MC simulation model was created for the Mediso AnyScan SPECT system based on specifications provided by the manufacturer and physical measurements. GATE version 8.2 with Geant4 version 10.05 was used. A detector head model was created and positioned for either dual or triple head configuration, as per the physical camera. Each head contained a 9.5 mm NaI thick crystal set as a ‘sensitive volume’, meaning all interactions within the volume were recorded. Geometric structures for the glass light-guide and the compartment housing the signal electronics were modelled in the simulation, as they were found to contribute to scatter into the crystal (as seen in other work [23]). The individual photomultiplier tubes were not defined in the simulation; instead, a general ‘back-compartment’ material was used to mimic the density of the electrical components, as suggested in [3]. In addition to SPECT, the AnyScan SCP system has a computed tomography (CT) and a positron emission tomography (PET) ring, but these were not included in the simulation.

Models for LEHR and MLEGP collimators were defined as uniform blocks of lead with a rectangular array of hexagonal holes according to technical specifications of hole size and septal thickness, shown in Table 1. An aluminium touch-plate was included in front of the collimator and the detector head was surrounded by lead shielding. The plastic casing surrounding the heads was not included in the simulation. A visualisation of the simulation model with a geometric cylindrical phantom on the patient bed is shown in Fig. 1, along with a photograph of the physical scanner.

The patient bed was defined from a CT image, in order to accurately describe its internal structure.

The CT scan of the phantoms was used to define the source geometry and attenuation distribution in the simulation. The Hounsfield unit (HU) value of each voxel was attributed to a material density and elemental composition. The CT images were manually segmented to determine regions for activity. The CT scan was used to ensure accurate replication of the positioning and material composition of the phantom in the simulation. The ‘ImageNestedParameterisedVolume’ feature in GATE was used to speed-up the simulation of voxelised volumes [24]. For the calibration ampoules and point sources, geometric source and phantom distribution were used due to the more simple geometry and material composition.

Nuclear data, physics and digitisation models

Each radioactive source was defined as an isotropic ‘UserSpectrum’ with nuclear data on radioactive emissions, energies and half-life from the ENSDF database [22]. All emissions of \(\upgamma\) rays and X-rays were included as a discrete energy histogram according to the energy and intensity provided by ENSDF. For \(^\)Lu, the \(\upbeta\) particle emissions were also included in the source definitions as an ‘Arb’ histogram with linear interpolation, so that any contributions from Bremsstrahlung were included. Radioactive decay was included in all simulated acquisitions. In each case, the full experimental activity was simulated.

All physics processes such as scatter and attenuation within the phantom, collimator and detector system were defined with the Geant4 emstandard_opt4 library with its default parameters; this physics library has undergone extensive validation [25]. Production range cuts of 0.05 mm were applied to photon tracking and 0.01 mm to electron tracking in the simulation. These are converted into an energy threshold for each material, below which no secondary particles are generated.

The scintillation process and light detection were not incorporated into the simulation due to large computation times. Instead, the photon-detection process was modelled through GATE digitisation modules, which convert photon interactions in the crystal into digital counts and apply blurring to mimic the response of the physical system. An ‘Adder’ module was used to sum detected events. A Gaussian blurring of the energy spectrum was applied by with the ‘blurring’ command, defining an energy resolution of 9.5% at 140 keV, as stated by the manufacturer. The energy resolution was extrapolated to other energies using GATE’s in-built default inverse square law function. A Gaussian intrinsic spatial blurring was defined with GATE’s ‘spblurring’ with full-width half maximum of 3.2 mm as quoted by the manufacturer [15]. The timing resolution of the detector was not modelled as the dead-time percentage was found to be insignificant at the activities used in this work (\(<0.05\%\) at 200 MBq of \(^}\)Tc which is almost twice the phantom activities used here).

Simulation execution and post-processing

A simulation was conducted for each of the experimental acquisitions of point-sources, calibration ampoules and phantoms. All simulations were split into parallel jobs and run on a HTCondor cluster [26] with 660 total CPU cores.

The GATE ‘Singles’ output of all simulations were recorded in TTrees in ROOT, an object-oriented programme developed by the CERN community [27]. Simulated energy spectra were extracted directly from the TTrees.

To assess the intrinsic efficiency, a Gaussian fit was applied to the photopeaks in the simulated energy spectra of the calibration ampoules. Equation 1 was used to calculate the absolute detection efficiency of the simulation as a function of energy; these are compared to the same measurements from the experimental data in Fig. 4a. The intrinsic efficiency measurements showed that the simulation consistently overestimated the efficiency and this effect was worse at low energy, suggesting that the simulation digitisation modules do not accurately represent the signal-processing chain of electronics in the physical system. To correct for this, the ratio of the experimental and simulated efficiency was calculated for each peak of the calibration sources. An efficiency correction function, \(\epsilon _}\), was modelled as a function of energy, E, using the equation

$$\begin \epsilon _}(E) = a + b \text(E) + c \text^2(E) \end$$

(3)

where a, b and c are the fitted parameters of the function and ln is the natural logarithm. Such a function is commonly used for the efficiency of NaI detectors [28]. The fitted \(\epsilon _}\) returned zero efficiency at 26 keV, rather than 2 keV seen in the experimental data. Thus, the efficiency was fixed to zero at 2 keV and a linear extrapolation was applied from the lowest-energy experimental data point. This energy region is significantly lower than energies typically used for clinical SPECT imaging, so the linear extrapolation was deemed sufficient. Figure 4b displays the logarithmic and linear functions applied to the efficiency-correction data. Residuals of \(\epsilon _}\) to the data above 46 keV (the lowest experimental data point) are shown in Fig. 4c; the error bars are the standard uncertainty on each data point (calculated with Eq. 2) and the dashed lines represent the uncertainty of the fit. The fit uncertainty was taken as the standard deviation of 1000 fits after perturbing each data point by a random number following a normal distribution with sigma equal to the point’s standard uncertainty.

A combination of the logarithmic and linear functions was used to correct the simulation efficiency:

$$\begin \epsilon _} = \left\ - 0.0134 + 0.0067 ~ E & \quad \text E < 46~}, \\ -3.755 + 1.702 ~\text(E) -0.168 ~\text^2(E) & \quad \text E \ge 46~\text . \end \right. \end$$

(4)

This function was applied to the simulation output through Poisson sampling [29] on an event-by-event basis to correct the efficiency of the simulation for all acquisitions. The experimental background energy spectra were time-corrected and added to the simulated energy spectra in ROOT. The manufacturer-provided software to extract the experimental energy spectra resulted in a slightly different number of counts compared to those from the projection images. The simulated energy spectra were therefore normalised to the same total counts as the experimental energy spectra for each phantom. This normalisation was applied only for the comparison of energy spectra and not in other analyses, so did not affect the sensitivity comparison.

To create projection images, ROOT was used to split each simulation into 120 projections to replicate each static position of a detector head and set gates for the relevant energy and scatter windows. An in-house code was then used to generate an interfile image for each projection, applying the same matrix and pixel sizes as the experimental acquisitions. The background projections with the relevant energy window were time-corrected (with Poisson sampling) and added to the simulated projections on a pixel-by-pixel basis.

Validation metrics

Several experimental observables were identified and used to validate the MC simulation. These were the sensitivity, energy spectrum, spatial resolution, sinograms (counts per time frame), and projection images.

The spatial resolution of the system was measured using the planar acquisitions of a point source filled with \(^}\)Tc in air. A linear profile with a 1-pixel width was defined through the centre of each simulated and experimental image, and the counts along that profile were recorded. The SciPy Orthogonal Distance Regression Python package [30] was used to apply Gaussian fits to each profile, accounting for uncertainty in both dimensions. The Poisson statistical uncertainty was used as the uncertainty on the counts for each pixel.

The other metrics were validated through tomographic phantom acquisitions of \(^}\)Tc. Manufacturer-provided software was used to extract the energy spectrum from all acquisitions for a direct comparison to the simulation. Sinograms were used to validate that the movement of the SPECT heads was replicated in the simulation and that the phantom and source had been defined appropriately. The peak signal-to-noise ratio was used as quantification of the similarity of simulated and experimental images. Peak signal-to-noise ratio, PSNR, is defined for comparison of a simulated image g to an experimental image f, both of size \(M \times N\), as

$$\begin \text (f,g) = 20 \log _ \left\_f } } } \right\} \end$$

(5)

where the mean squared error (MSE),

$$\begin \text = \frac \sum _^\sum _^ ( f(i,j) - g(i,j) )^2 \end$$

(6)

and \(\text _f\) is the maximum possible pixel value of image f [31]; all images were unsigned 16-bit, so this value was 32768. The PSNR value tends to infinity as the mean square error decreases, hence larger values of PSNR imply more-similar images. Since PSNR is not an absolute quantity, a reference was created to provide a baseline against which to compare the PSNR of the simulation and experiment. The reference PSNR was calculated by comparing the experimental images to themselves shifted by 1 pixel in each horizontal axis. The PSNR was calculated independently for each image slice, and the weighting mean and standard deviation of all slices was computed.

The absolute sensitivity, S, was calculated for the \(^}\)Tc and \(^\)Lu phantom acquisitions from the counts in each experimental and simulated image as

$$\begin S = \frac} \end$$

(7)

where c is the total counts in the image, A is the source activity, and T is the total acquisition time. This sensitivity was specific to the isotope, not the single emission as both emission and scatter windows were considered. Therefore, the branching ratio of specific photon emissions was not included in the absolute sensitivity calculation. Since the efficiency correction was applied to the simulation, its uncertainty was incorporated into the uncertainty of the sensitivity. The average uncertainty on the efficiency correction within the relevant energy window, \(\bar_}}\), was used. The uncertainty on the sensitivity, \(\sigma _S\), was calculated through

$$\begin \left( \frac\right) ^2 = \left( \frac\right) ^2 + \left( \frac\right) ^2 + \left( \frac_}}}_}}\right) ^2 \end$$

(8)

where \(\sigma _c\) and \(\sigma _A\) are the uncertainties on the counts and activity and the uncertainty on T is assumed to be negligible. For the simulation, the uncertainty on the activity was zero. The simulation uncertainty was dominated by the uncertainty on the efficiency-correction function; for the \(^}\)Tc and \(^\)Lu emission and scatter windows, \(\frac_}}}_}}\) ranged from 2.0 to 4.4%. For the experimental images, the efficiency term was zero and the uncertainty was dominated by the uncertainty in activity (which was below 1% for both isotopes).

The uncertainties on experimental and simulated data were quantified in each step of the validation process. Standard uncertainties (k=1) are quoted throughout this work unless otherwise specified.