SPECT/CT measurement of NEMA SPECT phantom to set up monte carlo simulations

A SPECT/CT measurement with a cold NEMA SPECT Phantom (i.e., the background compartment was filled with non-radioactive water) was performed at the University Hospital Würzburg to extract the most important phantom features for generation of a computerized NEMA Phantom. The measurements were performed with a Siemens Intevo Bold SPECT/CT system (9.5 mm crystal, medium-energy low penetration collimator). SPECT was performed to measure a realistic detector orbit (180° detector configuration, automatic contouring, continuous mode, 2 × 60 views). A high resolution CT (CT-HD, CTDIVol = 3.12 mGy, 28 mAs, 130 kVp, pitch factor 1.5, slice thickness 1 mm) was performed to determine the sphere positions inside the phantom. In addition, a standard CT (CT-AC, CTDIVol = 2.90 mGy, 28 mAs, 130 kVp, pitch factor 1.5, slice thickness 3 mm) was performed to generate a realistic mass-density map. Both CT images had an in-plane pixel size of 0.98 × 0.98 mm2 with a matrix size of 512 × 512 pixels. Based on the CT-HD, the sphere centers were determined using a semi-automatic thresholding method, the details of which are explained in the Supplementary Material.

Monte carlo simulations for all possible sphere permutations

Monte Carlo (MC) simulations were performed in SIMIND [22] to determine the influence of sphere positioning on the RCs. Figure 1 gives a schematic overview of the workflow used to generate the dataset used in this study. First, CT-AC was interpolated to a voxel size of 2.4 × 2.4 × 2.4 mm3 (matrix size of 256 × 256 × 256). Next, the HU values of the CT were converted to mass density by using a pre-determined, scanner-specific two-segment linear function [23]. The mass density map was then used as input for SIMIND to correctly simulate photon attenuation. The distance between detector and image center was extracted from the DICOM header and used for all performed simulations. The simulations were performed using SIMIND’s multiple sphere routine, where sphere centers \(_\) and diameters \(_\) are given as an input. The sphere diameters were permuted to create all 720 possible sphere configurations at the six designated positions, resulting in 720 simulations. In theory, mirroring the phantom in the sagittal plane along the gantry axis would represent an equivalent spherical arrangement, which would only require 360 simulations. However, the trajectory is not axisymmetric due to a rail system below the bed, which is why all 720 sphere permutations were performed. Examples of permutations of the computerized NEMA SPECT Phantom are given in Fig. 2. Strictly speaking, a separate mass density map would have to be created for each combination. However, the same mass density map can be used for all permutations without meaningful inaccuracies in the simulated attenuation, as the walls of the spheres are thin (around 1 mm according to manufacturer specification), and the attenuation coefficient is similar to that of the surrounding water (60 HU Vs.0 HU).

Fig. 1

Schematic overview of the generation of the dataset utilized in this study

Fig. 2

Example permutations of the computerized NEMA SPECT Phantom

The simulation was set to recreate a measurement on the SPECT/CT system described above (Siemens Intevo Bold) using a 20% main energy window at 208 keV and two adjacent 10% scatter windows. In addition to the NEMA SPECT Phantom, simulations of the NEMA PET Phantom were also performed by replacing the 60 mm sphere by the 10 mm sphere. This resulted in a total of \(2\times 720=\text\) simulations, which were performed on the local high-performance computing cluster at the University of Würzburg. The variance-reduction in the SIMIND program makes the noise from the Monte Carlo simulations not following a Poisson distribution. Hence, a large number of histories were run in order to achieve convergence in simulated projections. After scaling the projections (activity concentration 2 MBq/mL), Poisson noise was added to generate realistic projections.

SPECT reconstructions of simulations

In order to perform reconstructions for such a large number of simulations in a time-efficient manner, open-source reconstruction programs that provide batch reconstruction capabilities were used. Ordered Subset Expectation Maximization (OSEM) reconstructions without RR (OSEM_noRR) were performed using CASToR [24] (attenuation correction (AC), scatter correction (SC) with triple energy window (TEW) method, voxel size: 4.8 × 4.8 × 4.8 mm3, matrix size: 128 × 128 × 128) with 10 subsets and an increasing number of iterations (1 to 10 with step size 1, corresponding to 1 to 100 updates). Additionally, OSEM reconstructions with RR (OSEM_RR) were performed using STIR [25] (AC, SC with TEW, voxel size: 4.8 × 4.8 × 4.8 mm3, matrix size: 128 × 128 × 128) with 10 subsets and an increasing number of iterations (2 to 20 with step size 2, corresponding to 2 to 200 updates), a reconstruction very similar to the Siemens Flash3D reconstruction. Both reconstructions were performed for all \(2\times 720=\text\) permutations of the NEMA PET Phantom and the NEMA SPECT phantom, respectively.

Determination of recovery coefficients of the permutations

To determine a set of six RC values for each configuration, both SPECT reconstructions (OSEM_noRR and OSEM_RR) were interpolated to CT-HD resolution using tri-linear interpolation. Although linear interpolation can lead to a small degradation in recovery, it was applied here in order to compare recovery using the high-resolution activity mask. An alternative approach would have been to interpolate the segmentation mask to SPECT resolution. However, this would have resulted in a more pronounced sampling effect due to the larger SPECT voxels. Then, a set of RC values was calculated for both phantoms, both reconstruction types, and all numbers of iterations (10 different for OSEM_noRR and OSEM_RR), as described in the Supplement. The RC is defined as the SPECT-based total activity in the sphere divided by the nominal activity within the sphere. For each sphere the mean RC (\(\overline\)), the maximum RC (\(_\)) and the minimum RC (\(_\)) over all 720 permutations were determined.

To quantify a potential spread of RC for each sphere, the RC variation \(_\) was calculated as the difference between \(_\) and \(_\) divided by \(\overline\) (Table 1):

Table 1 Statistical analysis of the RCs for all 720 permutationsFit of the RC curve

For each of the 1,440 simulations, a non-linear least squares fit of the RC curve (RC values plotted against the sphere diameter) was performed for both SPECT reconstructions (OSEM_noRR, OSEM_RR) using the curve_fit function ( [26], non-linear least squares fit) of the Python package scipy, with the following fit function [27, 28]:

$$_\left(d\right)=\right)}^\right)}^$$

(2)

For each fit, the values of \(\beta\) and \(\gamma\), as well as the coefficient of determination \(^\) were determined. A mean fit \(\overline_}\left(d\right)\) was then determined over all 720 permutations of each phantom type as follows:

$$\overline_}\left(d\right)=\frac\sum _^__}\left(d\right)$$

(3)

The standard deviation \(_\left(d\right)\) was calculated using the following formula:

$$_\left(d\right)=\sqrt_^__}\left(d\right)-\overline_}\left(d\right)\right)}^}$$

(4)

Just as for RC (see Eq. 1), the variation of \(\beta\) or \(\gamma\) over all permutations is determined by the following formulas:

$$_=\frac_-_}}\,\,\, \text\, \,\,_=\frac_-_}}$$

(5)

In order to determine the global goodness of the fits to the determined RC values, the coefficient of determination \(^\) was calculated for each fit performed.

SPECT measurements for different sphere positioning

In order to test the simulation-based observations under real measurement conditions, additional 177Lu SPECT/CT measurements were performed on a NEMA SPECT Phantom (different copy of the same phantom) at the University Hospitals Leuven. These measurements were based on the preliminary accreditation protocol given by EARL. In total, six different sphere configurations were measured. A Standard set of sphere configurations, where the spheres were positioned with alternating sphere sizes (large sphere surrounded by the two smallest spheres; this positioning is proposed in the preliminary EARL accreditation protocol (cite)) and then the entire sphere insert was rotated by 0°, 120°, and 210° (Standard_0, Standard_120, Standard_210); and a second set of sphere configurations (Switch) in which the smallest to the largest sphere were arranged clockwise and again rotating the entire sphere insert by 0°, 120° and 210° (Switch_0, Switch_120, Switch_210). An overview of all six sphere arrangements is given in Fig. 3. The measurements were performed with the same camera model as the one reproduced in the simulations (Siemens Intevo Bold with 9.5-mm crystal thickness and medium-energy low-penetration collimator) using the same imaging (180° detector configuration, 60 projections for each detector, time per projection of 20 s (this resulted in the 25 kcts per detector and projection proposed in the preliminary EARL protocol), matrix size: 256 × 256, pixel size: 2.4 × 2.4 mm2, a 20% main energy window at 208 keV and two adjacent 10% scatter windows) and reconstruction parameters (Flash3D with RR, AC, SC with TEW, voxel size: 2.4 × 2.4 × 2.4 mm3, matrix size: 256 × 256 × 128, 2 subsets, and an increasing number of 5, 10, 15, 20, 25, 30, 40 and 80 iterations). From a batch of no carrier added lutetium trichloride (177LuCl3, EndolucinBeta, Isotope Technologies Munich), a syringe containing a target activity of 400 MBq was prepared using a well-type activity meter (–model VIK-202 with IBC-LITE software, Comecer). This activity was used to prepare a 200 mL stock solution from which the spheres were filled, resulting in a uniform activity concentration of 2 MBq/mL. To prevent the lutetium from sticking to the sphere walls, approximately 1 g of ethylenediaminetetraacetic acid was added to the stock solution. A high-resolution CT (CT-HD, CTDIVol = 10.20 mGy, 150 mAs, 110 kVp, pitch factor of 0.8, axial resolution of 1 mm) was acquired to determine the sphere positions, and a low-resolution CT (CT-AC, CTDIVol = 10.20 mGy, 150 mAs, 110 kVp, pitch factor of 0.8, slice thickness of 5 mm) was acquired for attenuation correction. Both CT images had an in-plane resolution of 0.98 × 0.98 mm2 with a matrix size of 512 × 512 pixels. The segmentation masks used to calculate the RCs were obtained using the same methodology as applied for the creation of the activity masks for the simulations (NEMA SPECT Phantom measurement at the University Hospital Würzburg). The RCs for all six sphere configurations were calculated for all different numbers of Flash3D iterations using Eq. 3 of the Supplementary Material. In addition, a fit of the RC curve was performed for each of the six measurements using the fit function given by Eq. 2.

Fig. 3

Overview of the six different sphere configurations used for the SPECT/CT measurements at the University Hospitals Leuven. The figure depicts the axial slice of the SPECT/CT fusion based on CT-HD and Flash3D SPECT reconstruction for all measured sphere configurations. Top from left to right: Standard_0, Standard_120, and Standard_210. Bottom from left to right: Switch_0, Switch_120, and Switch_210

Comparison between measurements and simulations

The agreement between the effects of sphere positioning as predicted by MC simulations and physical-phantom measurements was investigated. For this purpose, the RC curves of the six measured sphere configurations and of the simulations with the same sphere configuration were calculated. Due to an additional 30° rotation of the entire sphere insert against the sphere positions used for the simulations, the sphere configurations Standard_210 and Switch_210 were not present in the 720 simulated sphere permutations of our analysis. Therefore, additional SIMIND simulations were performed for these two configurations using the same parameters as the previous simulations. Given the visual similarity of the simulated and measured curves, it can be reasonably assumed that predictions about the specific RCs for each sphere configuration can be made based on the simulations. Additionally, two metrics were calculated to assess the deviation between measurements and simulations quantitatively:

$$\left(\text\right)= \frac_^\left|R_\left(d\right)-R_\left(d\right)\right|$$

(6)

$$\left(\text\right)= \frac_^\left|R_\left(d\right)-}_\left(d\right)\right|,$$

(7)

where \(j\) is one of the six measured configurations, \(R_\left(d\right)\) is the RC curve of the measurements and \(R_\left(d\right)\) is the RC curve of the simulations of these configurations. \(}_\left(d\right)\) is the mean RC curve over all six measured configurations.