PET image reconstruction using physical and mathematical modelling for time of flight PET-MR scanners in the STIR library

This work demonstrates how computational and physical modelling of the positron emission tomography (PET) image acquisition process for a state-of-the-art integrated PET and magnetic resonance imaging (PET-MR) system can produce images comparable to the manufacturer. The GE SIGNA PET/MR scanner is manufactured by General Electric and has time-of-flight (TOF) capabilities of about 390 ps. All software development took place in the Software for Tomographic Image Reconstruction (STIR: http://stir.sf.net) library, which is a widely used open source software to reconstruct data as exported from emission tomography scanners. The new software developments will be integrated into STIR, providing the opportunity for researchers worldwide to establish and expand their image reconstruction methods. Furthermore, this work is of particular significance as it provides the first validation of TOF PET image reconstruction for real scanner datasets using the STIR library. This paper presents the methodology, analysis, and critical issues encountered in implementing an independent reconstruction software package. Acquired PET data were processed via several appropriate algorithms which are necessary to produce an accurate and precise quantitative image. This included mathematical, physical and anatomical modelling of the patient and simulation of various aspects of the acquisition. These included modelling of random coincidences using ‘singles’ rates per crystals, detector efficiencies and geometric effects. Attenuation effects were calculated by using the STIR’s attenuation correction model. Modelling all these effects within the system matrix allowed the reconstruction of PET images which demonstrates the metabolic uptake of the administered radio-pharmaceutical. These implementations were validated using measured phantom and clinical datasets. The developments are tested using the ordered subset expectation maximisation (OSEM) and the more recently proposed kernelised expectation maximisation (KEM) algorithm which incorporates anatomical information from MR images into PET reconstruction.


Introduction
Positron emission tomography (PET) is an important tomographic imaging modality with particular deployment in cardiology [1], neurology [2] and oncology. The successful utility of PET in these branches of medicine can be seen by its increasing use for various clinical indications. PET scans are non-invasive and helpful in early diagnosis and management of conditions such as in Parkinson's disease [3], dementia [4,5], coronary artery disease [6] and oncology [7] by quantitatively assessing metabolic and functional alterations [8]. PET scans are also clinically useful in monitoring the progression of diseases and the patient's response to therapy [9].
Being of interest in both the clinical and research communities, a platform that can readily provide an open framework to reconstruct https://doi.org/10.1016/j.ymeth.2020.01.005 data is of high significance. The reconstruction software used in this work is an open source library based on object-oriented programming called 'Software for Tomographic Image Reconstruction' (STIR: http:// stir.sf.net) [10]. STIR is a widely used C++ library providing tools to reconstruct emission tomography images with all data corrections such as normalisation, attenuation, randoms and scatter estimation. STIR also allows manipulation of projection data and images, image analysis and analytical simulations. It is designed to be modular and has an inheritance between classes. This allows class extension and thus can incorporate novel algorithms and scanner geometries conveniently without the need to re-implement existing attributes and functions. This structure provides a platform to the PET research community where the library is expandable to "fit the needs" of new users as the field progresses by deriving new classes from basic building blocks. For instance, to implement novel iterative reconstruction algorithms, an existing reconstruction base class such as Ordered Subset Maximum A Posteriori One Step Late can be extended and other STIR functionalities can be reused [11]. Recent collaborative efforts within the STIR developer's community have led to the extension of the library to novel iterative reconstruction algorithms such as the kernel expectation maximization (KEM) [12] or modelling timing information for time-of-flight (TOF) PET [13].
PET scanners measure data which are difficult to reconstruct, especially with sufficient accuracy without knowing all necessary information provided by the manufacturer. Despite the wide range of novel iterative reconstruction algorithms available within STIR, it can be challenging for users to accurately reconstruct data measured from clinical or preclinical PET scanners. One needs to carefully implement proprietary scanner data formats accurately within STIR as it currently supports only CTI ECAT6, ECAT7 and Siemens Biograph mMR file formats.

Objectives
The aims of this work are: 1. Modelling of the GE SIGNA PET/MR scanner's acquisition process. 2. Expand STIR library to GE formats. 3. Validation of the images reconstructed by STIR and the vendor's software. 4. Demonstration of reconstruction using KEM. 5. Demonstration of TOF reconstruction using STIR library for measured data.

GE SIGNA PET/MR scanner in STIR
GE SIGNA TOF PET/MR scanner installed at Invicro, Imperial College London, Hammersmith Hospital, London, UK, was used in this study. It is a 3T MR scanner that has an integrated MR-compatible PET scanner [14]. The PET scanner has a field of view of 60 cm and 25 cm in transaxial and axial directions, respectively. It is composed of lutetiumyttrium oxyorthosilicate (LYSO) crystals, silicon photomultipliers (SiPM) and other processing electronics [15]. The SiPM together with the LYSO crystals allow for TOF imaging [16]. The PET component of GE SIGNA comprises of 28 detector modules and each detector module has 5 detector units (each detector units are separated from the next unit by a crystal-free gap of size 2.8 mm). Further, there are 4 blocks per unit and each block has × 4 9 crystals. Each block further has a × 3 6 array of SiPMs. The scanner is configured by detector modules placed next to each other forming a cylindrical bore with detector units along the axial direction. Thus, the scanner has 45 detector crystals along the axial direction and 448 detectors around the scanner making a total of 20,160 detectors. The modules are placed with module 0 centered above the isocenter of the scanner which implies that crystal 0 is not vertically above the center of the scanner. Crystal 0 is located at an offset of°5.23 from the vertical center when viewing from the back of the scanner. The LYSO crystals are 5.3 mm along the axial direction, 3.95 mm along the transaxial direction and 25 mm deep. With the fast detection electronics of the scanner, its timing resolution is 390 ps. The scanner has a coincidence timing window of 4.57 ns and timing least significant bit of 13.02 ps duration which leads to 351 timing bins. The scanner saves these as signed timing bins (i.e. numbered from −175 to + 175 bins) which the software combines in groups of 13 timing bins, producing 27 sinogram timing bins. The coincidence events are stored in a sinogram of dimensions × × × 224 1981 27 357, where 224 corresponds to the total number of views, 1981 corresponds to the total number of axial positions, 27 corresponds to the number of "mashed" timing bins and 357 corresponds to the number of tangential positions. The number of views refers to half of the total number of detectors along the ring by convention.
The scanner parameters are included in the scanner template class "Scanner.cxx" [17]. This class stores the relevant geometric information for PET scanners.

Scanner model
The radioactive distribution of the object placed within the scanner is mathematically modelled as follows: where y represents the acquired sinogram measurement, A represents the system matrix where x represents the image of the radioactivity distribution of the object and b represents the background coincidences (e.g. randoms and scatter). The system matrix models the sensitivity at each voxel m of the scanner by modelling the probability of positron annihilated in voxel m and being detected as a photon-pair coincidence by the detector pair p. The sensitivity of each voxel m combines the detector efficiency, geometric and attenuation effects.
Specifically, the TOF information for the measured sinogram data is modelled within the forward model as: (2) where y pt represents the time of flight projection data measured by the scanner and the system matrix a pt m ; models the probability of event emitted in voxel m to be detected by detector pair p within the signed timing bin number, t. b pt corresponds to the background counts of the corresponding timing bin t and detector pair p.
In this paper, the well known ordered subset expectation maximisation (OSEM) [18] iterative image reconstruction algorithm is extended to include timing information (i.e. TOF-OSEM) and MR information (i.e. TOF-KEM). These algorithms are used to reconstruct the measured sinogram y pt as measured by the scanner. Currently only sinogram-based reconstructions are demonstrated.

TOF-OSEM algorithm
The TOF-OSEM algorithm was described in [13]. Although the TOF expansion of the STIR library was recently implemented and validated on simulated datasets [13], this is the first time that the implementation is tested on acquired datasets.
The TOF-OSEM algorithm is described as follows: , st represents the correction term to account for the sum of randoms and scatter events (or as will be referred to background term) for subset s P. Wadhwa, et al. Methods xxx (xxxx) xxx-xxx and timing bin t and S j represents the subset of detector pairs. Note that one full iteration, or simply mentioned as iteration, corresponds to a full cycle of sub-iterations across all subsets (i.e. iterations = sub-iterations/subsets).

TOF-KEM algorithm
The KEM algorithm was first proposed by Wang et al. [19] and was developed in STIR library by Deidda et al. [12,20]. This has been now extended further to accommodate TOF data as well. The algorithm assumes that an image is described by the sum of a number of weighted kernels, and it tries to reconstruct the weight of each kernel by assuming it knows the kernel from another source of information (e.g. MRI): where value at voxel m of image x m represents the linear combination of the kernel matrix (k fm represents its fm th elements), N m represents the number of feature vectors and kernel weights . The kernel weights are reconstructed as follows: where k lm represents the MR kernel image. The final image which contains the information of interest (i.e. the radioactivity distribution) is obtained by applying Eq. (4) after the last iteration of Eq. (5).

GE SIGNA PET/MR scanner acquisition data in STIR
Image reconstruction of data acquired with the GE scanner using STIR library requires the incorporation of input and output information that convert the data into an accessible data format. Specifically, the output files saved by the scanner console at the end of a PET scan include: 1. Listmode (LM) and sinogram raw data files (RDF) compressed with the GE proprietary compression algorithm 2. Raw data calibration files including normalisation, geometric and well counter calibration factors 3. PET Images For Attenuation (PIFA) files (i.e. processed PET attenuation correction (AC) images obtained by converting MR images acquired using MR sequences including zero echo time and Dixon) 4. Reconstructed PET images with reconstruction settings selected at the time of scan All the output raw data files are in GE proprietary data format and can be processed offline solely with the GE's proprietary reconstruction software. GE stores data and header information in Hierarchical Data Format 5 (HDF5) file format.
HDF5 files include: 1. Sinogram RDF generally called 'rdf.1.1' 2. LM file generally called 'LIST0000.BLF' 3. Normalisation data file, 'norm3d', and geometric correction data file, 'geo3d' These files are used as input to validate the implementation of the acquisition model for this scanner in STIR.
However, STIR has its own native sinogram data format, used as input for the iterative and analytical reconstruction algorithms. Thus, all the HDF5 files need to be converted into STIR native interfile format [17]. A wrapper class was created to import the HDF5 files during this work. HDF5wrapper class has functions implemented to access and read the required information which is further modelled within the system matrix as discussed in Section 2.4.

TOF Acquisition Data in STIR
TOF information is implemented in STIR within the system matrix as: This equation models the TOF information as a time spread function, T ptm for each time bin t and separates it from the sensitivity information, a pm . The implementation has been recently presented and validated in simulated data [13], but not yet for acquired data.
The scanner saves the detector pair and timing information in LM and sinogram RDFs. The data saved in the LM HDF5 file is accessible by navigating to the field called '/ListData/ListData' [21,17]. LM files store each event of an acquisition as 6 byte structure. This includes information such as the crystal identifiers (IDs) which registered the coincident photons and the timing bin. The detector pair positions are handled as axial and transaxial coordinates within the "CListModeDa-taGESIGNA" class of STIR. The signed TOF bin is handled within the "CListRecordGESIGNA" class. This is done by reading the TOF information along with the TOF "mashing" information and use TOF projection data, y t within the reconstruction framework. In the current implementation, TOF bins can be mashed together to get non-TOF data. The data is rearranged to be accurately read in STIR by applying transformations from scanner space to STIR space.
The next step to reconstruct PET images with STIR requires the calculation of the sensitivity image. This includes modelling of detector efficiency, geometric and attenuation factors. The sensitivity images demonstrate the spatial variation of expected coincidences per unit activity concentration.

Implementing normalisation correction for GE SIGNA in STIR
The geometric and detector efficiency factors are multiplied to get the normalisation correction factors per detector bin as below [22,23]: where n XY represents the normalisation correction factor for the detector pair, j which is comprised of detector IDs X and Y , X and Y represents the crystal efficiency factors, and g XY represents the geometric correction factors. The scanner exports HDF5 files, which are then used as inputs for calculating the normalisation factors [17]. Crystal efficiency factors are stored in 'norm3d' and geometric factors are stored in 'geo3d'. The norm3d HDF5 file stores crystal efficiency factors which are read as an array of × 448 45 in STIR. The 3D geometric correction factors are stored as a part of the projection data and only 16 views are stored i.e. for one module. These view datasets are of dimensions × 1981 357. The exported geometric factors are repeated 14 times within STIR to get full geometric correction projection data. Dead-time correction has not been implemented in this study.

PIFA-based attenuation correction
Attenuation correction was carried out using scanner extracted PIFA images. This was achieved by creating custom scripts that imported the PIFA images within the STIR reconstruction framework [17].
Furthermore, a background term needs to be calculated to correct for scatter and random events.
where R XY is the number of random events (or counts) detected by the detector pair X and Y , represents the coincidence time window (4.57 ns for this scanner), T represents the total acquisition time, S X and S Y represent the single events detected by crystal X and Y, respectively. Single rates for each crystal are exported from the LM HDF5 file and converted into random rates using the STIR utility 'con-struct_randoms_from_GEsingles' implemented during this work. Singles per second are detected for every crystal and then stored as a dataset with dimensions × 45 448 in LM HDF5 file. Each cell of the dataset corresponds to the number of singles detected by the crystal number X Y ( , ) where X corresponds to the ring number and Y corresponds to the transaxial crystal number. For an acquisition of time of N seconds, there is a list of + N 1 samples. There are + N 1 samples as the periodic data are asynchronous with the start of PET scan which needs to be taken into account during data calculation by weighting the first and last samples of the periodic data appropriately. Eq. (8) calculates the randoms correction sinogram, which is later used for every TOF bin after dividing it with the total number of timing bins. In the current STIR implementation, dead-time and isotope decay modelling are not considered within randoms correction.

Scatter correction
Non-TOF scatter correction sinogram was extracted from the vendor's reconstruction toolbox for scatter correction.
TOF-scatter is currently not implemented, whereas non-TOF scatter correction sinograms were approximated to TOF projection data with appropriate scaling factors for each TOF bins. The scaling factors were extracted using the emission projection data by calculating the fraction of counts in each TOF bin with respect to total number of counts.
The randoms and scatter correction TOF projection data were added to calculate the background term, b t .

Data acquisition, image reconstruction and analysis 2.5.1. Phantom acquisition
The following two phantoms were scanned: 1. Volumetric Quality Control (VQC) Phantom: Comprised of five 0.7 MBq 68 Ge spheres of 19 mm in diameter with 14 mm active diameter (with half mm tolerance) encased within MR-visible tubes made of NiCl 2 . The entire configuration is encased within a foam cube. The phantom was scanned for 10 min amd × 5 10 6 prompts were registered. 2. Hoffman Phantom: Scanned for 20 min and × 1.5 10 8 prompts were registered [26].

Clinical acquisition
A patient with radiotherapy induced pulmonary fibrosis [27,28], injected with 40.62 MBq of 18 F experimental radiotracer was scanned. The radiotracer was injected 90 min prior to the scan, which had a duration of 13 min. × 3.8 10 7 prompts were registered. MRAC images, PIFA images, uncompressed PET LM and RDF files were exported from the scanner for all above phantom and clinical acquisitions.

Image reconstruction
The scanner acquired uncompressed LM data which were histogrammed into TOF STIR sinograms using the existing 'lm_to_projdata' utility. This is possible due to the implementation described in 2.4.1.
Data correction sinograms were calculated using custom utilities and parameter classes with STIR for this scanner. Normalisation corrections were calculated using the 'correct_projdata' utility with parameters set to read data from 'norm3d' file exported from the scanner. Attenuation correction was carried out using the exported PIFA images as demonstrated in Fig. 1. Randoms correction sinogram was calculated using 'construct_randoms_from_GEsingles' utility. Scatter correction was done using the sinogram extracted from the vendor's reconstruction toolbox. For TOF image reconstructions, the TOF projection data were histogrammed using the LM file. Randoms and scatter correction sinograms were translated into TOF projection data by developing custom utilities and appropriate scaling of data as described above in sections 2.4.4 and 2.4.5. Attenuation and normalisation sinograms were uniform in timing bins and were multiplied with the background term to create additive sinograms for STIR implementations. Image reconstructions were carried out using OSEM and TOF-OSEM algorithms with and without point spread function (PSF) modelling. TOF-KEM reconstructions were also carried out for the patient dataset. The MR image shown in Fig. 1 was used to derive the kernels. All reconstructions were carried out with 28 subsets.

Image reconstruction with vendor's reconstruction toolbox
GE proprietary reconstruction algorithms including VUE Point HD (or fully 3D iterative reconstruction which is further referred to OSEM-GE), VUE Point FX (or fully 3D TOF iterative reconstruction which is further referred to TOF-OSEM-GE), VUE point HD SharpIR (or PSF-OSEM-GE) and VUE point FX SharpIR (or PSF-TOF-OSEM-GE) were used to reconstruct VQC, Hoffman and patient datasets with 28 subsets [29,30].

Image analysis
Emission and data correction sinograms exported from the scanner are calculated using the implementations made in the STIR and then compared using voxelwise subtraction. TOF-OSEM reconstructions with STIR and the vendor's reconstruction toolbox (or 'GE' as will be now referred to in this document) for the VQC phantom were used to compare the full-width-half-maximum (FWHM) [31]. FWHM was calculated using a STIR utility called 'find_fwhm_in_image', for all five spheres for vendor's reconstruction toolbox and STIR based reconstruction. The FWHM for all spheres in all three dimensions were averaged and standard error was calculated as the tolerance window. TOF-OSEM reconstructions with STIR and vendor's reconstruction toolbox for the Hoffman phantom and the patient dataset were compared using standardized uptake value ratio (SUVR) and coefficient of variation (CoV). SUVR is calculated as the ratio of the SUV of target and reference regions. In this study, a region of interest (ROI) drawn within the liver and spleen were used as the target and another ROI of same volume within the lung was used as the reference. CoV is calculated as the % of standard deviation over the mean. STIR reconstructions were also compared with the vendor's reconstructions using structural similarity index (SSIM) as described in [32]. The global SSIM values are reported as the similarity measure between STIR and GE reconstructed images and are calculated as:  Table 1 demonstrates the average of FWHM for all five 68 Ge spheres in all dimensions for the VQC phantom. The comparisons are demonstrated for images reconstructed with OSEM, PSF-OSEM, TOF-OSEM and PSF-TOF-OSEM algorithms with vendor's reconstruction toolbox and STIR over 28 subsets and 3 subsequent iterations. Tables 2 and 3 demonstrate the SUVR comparisons for liver and spleen ROIs respectively, compared with lung ROI. The comparisons are made for image reconstructed with patient dataset. Iterative algorithms compared are OSEM, PSF-OSEM, TOF-OSEM, PSF-TOF-OSEM, TOF-KEM and PSF-TOF-KEM. Fig. 2 compares the direct ring sinograms from STIR and vendor's reconstruction toolbox for segment 0. The comparison is made using the uncompressed listmode file for all the various phantom and clinical datasets but only demonstrated for VQC phantom here. The unlisted sinogram using STIR utility and vendor's reconstruction toolbox is demonstrated here for TOF bin 0 and 2, segment 0 and axial position 18. Fig. 3 compares the non-TOF sinograms from STIR and vendor's reconstruction toolbox. The comparisons were made for emission, normalisation and randoms sinograms. The relative difference has been demonstrated by performing voxelwise subtractions. Fig. 4 and 5 demonstrates the non-TOF reconstruction with OSEM and PSF-OSEM respectively, for 28 subsets and 3 iterations for Hoffman phantom. Fig. 6 and 7 demonstrates the non-TOF reconstruction with OSEM and PSF-OSEM respectively, for 28 subsets and 3 iterations for clinical dataset. The demonstration compares the reconstructions between STIR and vendor's reconstruction software visually.

Discussion
The TOF emission projection data extracted with STIR for all phantom and clinical datasets in the described implementations is identical to the acquired emission data. The non-TOF emission sinograms extracted from STIR and vendor's reconstruction toolbox for the VQC phantom dataset, as can be seen in Fig. 3, are identical, whereas there are small differences for the normalisation and randoms sinograms of STIR implementation and compared with the vendor's reconstruction toolbox. The relative difference in the normalisation correction sinogram is ± × 1.6 10 7 . The relative difference for the randoms sinogram is 0.6% to 1.1%. The differences between the randoms and normalisation sinograms are due to the fact that no dead-time modelling has been implemented in the current STIR version. All the data acquisitions that are used in this study have relatively low activities as compared to clinical situations. This implies that the dead-time correction is very small. Also, the randoms correction implemented in STIR is not identical to that in the vendor's reconstruction toolbox calculation [24,33].
FWHM comparisons for the VQC phantom dataset, as demonstrated in Table 1, shows that the STIR reconstructions have similar resolution with the vendor's reconstruction toolbox. There is a relative difference of 5% to 13% in the measurements with STIR as compared to vendor's reconstruction software.
Hoffman and patient datasets reconstructed with the OSEM algorithm with STIR and the vendor's reconstruction software, as shown in Fig. 4 and 6, shows visual differences which are a consequence of absent detector gap modelling within the system matrix in STIR and a lack  Methods xxx (xxxx) xxx-xxx of deadtime modelling in randoms. The detector gaps are not taken into account as STIR currently supports only the cylindrical scanner geometries and without any gaps. There are also minor differences that arise from the offsets between the PET and MR gantry. GE converts the PET reconstructions in MR space which is not currently automated within STIR. Although, these offsets are manually taken into account during this work. SUVR comparisons for OSEM-STIR and PSF-OSEM-STIR for patient datasets are in good agreement with the vendor's reconstruction software counterparts for liver/lung comparisons, as can be seen in Table 2. For the liver/lung, there is a maximum of 3% relative difference in the SUVR values. For spleen/lung comparisons there is a greater relative difference of 21% for OSEM due to differences in projectors and random modelling. Similar patterns are observed in the Hoffman phantom dataset. The patient images reconstructed with the two software packages with TOF-OSEM algorithms appear to have differences due to the reasons aforementioned. TOF images had differences in noise which is due to the lack of implementation of TOF scatter simulation in the current version of the STIR library. SUVR comparisons for TOF reconstructions further shows the effect of lack of TOF scatter implementation as there are greater differences between STIR and vendor's reconstruction software quantifications. TOF-KEM SUVR comparisons display a very slight improvement in terms of the quantification.
It can be seen in Fig. 10 that PSF-TOF-KEM has the best uniformity for spleen and TOF-KEM has the best uniformity for liver. This shows the effect of the MR anatomical information in PET reconstructions. TOF-KEM CoV values improves as the number of iterations increments.
The global SSIM value is between 0.84-0.85 for all comparisons performed on reconstructed images over 6 iterations. This shows a high structural similarity between STIR and GE reconstructed images. The difference between two images is majorly due to intensity difference between the two images. This study does not take into account TOF scatter correction which gives rise to differences in noise properties between STIR and GE reconstructions. The randoms correction modelling implemented in STIR is not identical to that in GE toolbox. GE reconstructions are scaled by a global and a calibration factor that is not currently taken into account within the STIR reconstructions.

Conclusion
In this study, the acquisition process for GE SIGNA PET/MR scanner is modelled within the open source library, STIR. All validations were made by comparing images reconstructed with iterative reconstruction algorithms available in the library for acquired phantom and clinical datasets. Reconstructed images appear to be at a comparable level with the ones provided by the manufacturer. Dead-time correction, absence of detector gap modelling and TOF scatter correction have not been modelled in this work, and may account for the observed discrepancies. The main advantage of modelling the acquisition processes of this scanner in the STIR library is that the scientific community can perform their research in image reconstruction with complete control of how the data is manipulated. This will help to pave the way for harmonisation of methodology across various PET scanners used in the clinic. Future work will also include calculating scatter correction with STIR using single Compton scatter simulation [34] and to integrate it with the next release of SIRF (Synergistic Image Reconstruction Framework) [35].