PET rapid image reconstruction challenge (PETRIC)

1 Introduction

Positron emission tomography (PET) is an advanced imaging technique which uses photon pairs (produced by positron annihilation) to detect the approximate location of the positron emissions. Even with the current scanners’ fast PET detectors and electronics, the measured photon pairs still pose an inverse problem [1]. For this purpose, several tomographic image reconstruction algorithms have been developed over the years, first by Hounsfield and Cormack, joint recipients of the 1979 Nobel Prize for Physiology or Medicine. There are three (not necessarily distinct) categories of reconstruction algorithms [2]: analytical [3, 4], iterative (or statistical) [5], and artificial intelligence (AI) [6]. Iterative algorithms are currently preferred in clinical PET systems due to their detailed modelling of the acquisition processes (though AI-based algorithms are gradually being made clinically available). However, iterative algorithms can be computationally demanding. This is exacerbated by recent developments such as long axial field of view PET scanners, which provide much more data. Furthermore, some PET applications (such as dynamic imaging or motion correction) require subdividing acquisition data [7]. This results in multiple low-count reconstructions. The problem setting is therefore becoming increasingly challenging. In real practice, computational speedup in image reconstruction can permit running more complex (previously prohibitively slow) algorithms to provide significantly improved medical images.

Several advanced image reconstruction methods have been developed exploiting improved optimisation methods, machine learning, and more powerful computer hardware. However, it is difficult to objectively compare and assess such reconstruction algorithms and frameworks for accuracy and computational time.

Over the last two decades, several scientific communities have established the concept of “challenges” to facilitate fair comparison of different algorithms. These are well-known in the fields of image registration, segmentation, and reconstruction. Each challenge can focus on a particular topic in order to accelerate research in a specific direction. Inspired by these developments, in 2024 we held the PET Rapid Image Reconstruction Challenge: PETRIC. In this paper, we describe the design and rationale of the challenge.

2 PETRIC description

The primary aim of PETRIC was to stimulate research & development of faster PET image reconstruction algorithms applicable to real-world data. The challenge used a predetermined maximum a-posteriori (MAP, also known as penalised likelihood) problem and a “gold-standard” reconstruction algorithm which is known to converge to the optimum. Participants were given access to a sizeable set of phantom data acquired on a range of clinical scanners. The main task was to reach a solution close to the “gold standard” images (similar—within a threshold—mean radioactivity values within specified regions of interest) as quickly as possible (minimum computation time). One key component of the challenge was the use of the Synergistic Image Reconstruction Framework (SIRF) as a common computational platform for all participating reconstruction algorithms [8]. Additionally, the Core Imaging Library (CIL) was used to subclass its iterative algorithm class [9]. All algorithms therefore use the same projectors and objectives, and thus their performance is determined purely by algorithmic approach.

Seven teams registered for the challenge, with four teams submitting their work for evaluation. These four teams submitted a total of nine algorithms (each team contributing up to three entries).

The design, preparation and execution of PETRIC involved several coordinated major tasks:

• Gold standard definition
• Collection & curation of datasets
• Objective evaluation
• Definition of metrics & thresholds
• Live scoring & feedback
• Concluding workshop & awards

3 Materials and methods

3.1 Reconstruction problem

The optimisation problem is a MAP estimate using the smoothed relative difference prior (RDP) [10], i.e.,with y the measured data, x the reconstruction image andwhere the constraint set C is defined asfor some mask M provided. The parameter β>0 is the regularisation parameter that trades-off the fit to the data with the regularity provided by the prior. Note that SIRF defines the objective function as above, but algorithms from CIL minimise the function. The log-likelihood (up to terms independent of the image) iswith y a vector with the acquired data (histogrammed), and y^ the estimated data for a given image x,Here m denote the multiplicative factors corresponding to the detection efficiencies and the attenuation coefficient factors, a the “additive” background term corresponding to an estimate of randoms and scatter coincidences, pre-corrected with m), A is an approximation of the line integral operator [11], and diag an operator converting a vector to a diagonal matrix.

x*=argmaxx∈CΨ(x;y)

Ψ(x;y)=L(y;y^(x))−βR(x)(1)

x∈C⇔{xi≥0i∈Mxi=0else

L(y;y^)=∑kyklogy^k−y^k

y^(x)=diag(m)(Ax+a).(2)

Due to PET conventions, for some scanners, some data bins will always be zero, which corresponds to any potential “virtual crystals,” in which case any corresponding elements in m will also be zero. The corresponding term in the log-likelihood is defined as zero. All other elements in a are guaranteed to be (strictly) positive (ai>0).

The smoothed RDP is given bywith V the set of voxels, Ni the neighbourhood of voxel i (here taken as the 26 nearest neighbours), wi,j weight factors (here taken as “horizontal” voxel-size divided by Euclidean distance between the i and j voxels), κi an image to give voxel-dependent weights (here predetermined as the square root of minus the row-sum of the Hessian of the log-likelihood at an initial ordered subsets expectation maximisation (OSEM) [12] reconstruction, see [[13], Equation 25], γ an edge-preservation parameter (here taken as 2), and ε a small number to ensure smoothness (here predetermined from an initial OSEM reconstruction). The initial OSEM reconstruction was also given to all participants to potentially use as a starting point.

R(x)=12∑i∈V∑j∈Njwi,jκiκj(xi−xj)2xi+xj+γ|xi−xj|+ε,(3)

3.2 Reference reconstruction

Reference MAP solutions are computed using a modified version of the BSREM algorithm, using an EM-like preconditioner, subsets, and a relaxed step-size:with PC the projection onto the non-negative image space and S the average of the “sensitivity images” AsT⋅ms for each subset s, δ a small number (determined from the OSEM image) and α a step-size that decreased with iterations. This algorithm has been shown to converge to the MAP solution [14]. Our implementation [15] used a number of subsets different for each specific scanner, and manual tuning of the step size parameters. Nevertheless, we found that a high number of updates (of the order of 15 000) was often required. If necessary, the algorithm was restarted until our convergence metrics were satisfied, see Section 3.5. Two reference reconstructions for test datasets are shown in Figure 1. A dedicated example demonstrating PET image reconstruction with BSREM using SIRF can be found in a corresponding SyneRBI notebook.

x^new=PC(x+α(x+δ)S∇xsΨ(x;y)),

3.3 PETRIC datasets

3.3.1 Data collection

At the start of the challenge, we provided example acquired data of a small set of phantoms for participants to use for the development of their methods, complemented by further datasets during the challenge. Example test datasets are shown in Figure 1. An overview of the curated data is in Table 1. New datasets were acquired at different sites following the challenge training period, see Table 2. These datasets were used for scoring the performance of each algorithm. In this way the bias towards a certain vendor or scanner model was kept minimal. In addition, this approach ensured fairness of the competition because at the time of the challenge, nobody had access to the final ground truth data - including the organisers themselves.

Table 1. Training data made available to participants during the challenge.

Table 2. Evaluation data released after the submission deadline.

We requested datasets from several medical physics groups. All datasets had to be made open access, under the Creative Commons Attribution 4.0 International license. We encouraged the contributing scanning sites to provide not only the acquired raw data for the final testing but also the corresponding regions of interest. All phantom datasets have been made publicly available after the formal completion of the challenge. Participants were free to test their algorithms on additional datasets, for instance those available at the Zenodo SyneRBI Community. Human datasets were not included due to difficulties with organising data sharing agreements and privacy issues related to potential use of the data for facial recognition [16].

Only simple data from scanners and acquisition modes (e.g., only one bed position) were considered in this challenge. Reading the data was supported by the open source software for tomographic reconstruction [[17], STIR] either directly or through one of our conversion scripts. Some of these scripts have not been made publicly available yet due to pending approval from scanner vendors. The scanner models for which we can handle data are: All necessary datasets were saved in a STIR-supported version of the Interfile standard to allow offline reconstruction. A readme file with information about the phantom preparation and acquisition protocols was provided.

• Siemens Healthineers: Biograph mMR (native SIRF/STIR processing), Biograph Vision 600 (conversion via e7tools)
• GE HealthCare: Signa PET/MR (conversion via Duetto [18], but we have no access to the MRAC code), Discovery 690/710 (conversion via pettoolbox), Discovery MI (conversion via Duetto)
• Mediso: AnyScan (conversion handled by the National Physical Laboratory (NPL))
• Positrigo: NeuroLF from output in SAFIR file format [19]

The data collected included the following files: Curated datasets are available at https://petric.tomography.stfc.ac.uk/data/

• list-mode files, as exported from the scanner console,
• attenuation correction files as provided by the scanner (i.e., CTAC or MRAC) in DICOM, Nifti or STIR Interfile format,
• normalisation and calibration files, in the format supplied by the vendor,
• image reconstructed with vendor software (DICOM) with all corrections applied for cross-checking,
• volumes of interest (VOI) provided as images in exactly the same dimensions/voxel-size as the scanner reconstruction (Nifti or STIR Interfile format), one file for each VOI. VOIs were constructed by the sites contributing the data or constructed by the PETRIC team. Scripts were provided to construct VOIs for the NEMA Image Quality and Hoffman phantoms.

3.3.2 Data curation

To prepare the PET datasets for analysis, we followed a structured data curation workflow. The process was designed to ensure consistency, reproducibility, and compatibility with the SIRF framework. We describe this workflow here. More detail is provided in Appendix.

All datasets were organised under a standardised directory structure, where the dataset name followed the convention scannername_phantomname. Initial quality control checks were performed using a dedicated script to validate the integrity of the raw data. Subsequently, we generated initial images using one of the provided VOIs as a template. Since all VOIs shared the same geometry, any of them could be used interchangeably. In cases where VOIs were not available, we used either a vendor-supplied image or a placeholder. During image generation, some metadata—such as modality, radionuclide, and scan duration—was lost. We manually restored this information by editing the relevant image header files, referencing the original prompt data. For internal reference reconstructions, we also configured the number of subsets in the dataset settings file. To support automated figure generation, we specified the slices of interest in the PETRIC configuration files. The center-of-mass coordinates of the VOIs, obtained during quality control, were used to guide this selection. When VOIs were not pre-defined, we generated them ourselves, based on the OSEM images. After VOI creation, we re-ran the quality control process to produce additional plots and verify VOI alignment. To calibrate reconstruction parameters, we estimated the penalisation factor by comparing each dataset with a reference dataset acquired on the same scanner, followed by some manual tuning. For reference solutions, we performed BSREM reconstructions as mentioned above. Given the computational demands of this step, we monitored the reconstruction progress and image quality throughout. Once the reconstructions are completed, we adapted a plotting script to visualise BSREM metrics and copied the final reconstructed image to the designated reference location within the dataset directory. To finalise the dataset, we removed intermediate files and verified the presence and accuracy of documentation files such as README.md. The curated datasets were then transferred to a web server for distribution and further use.

3.4 Evaluation and scoring

We defined several metrics which were used to determine the fastest algorithm to reach target image quality, see next section. All algorithms were executed for each dataset until their solution reached the required thresholds on all selected metrics (for 10 subsequent updates) or their runtime exceeded one hour. For each metric, results from all teams were ranked according to the wall-clock time required to reach the threshold on our standard computational platform. Ranking is from worst (rank 1) to best (rank N), with best corresponding to “the fastest to reach threshold.” The overall rank for each algorithm was the average of all ranks for the individual metrics on each dataset. The algorithm with the highest overall rank won the challenge. Note that due to potential variability in the wall-clock timing as well as in using the stochastic algorithms, each reconstruction was executed three times, and the median of the wall-clock time was used for the ranking.

3.5 Metrics and thresholds

Each entry to the challenge was evaluated on metrics based on the object as a whole, the background and specific volumes of interest (VOI). They all use the voxel-wise root mean-squared-error (RMSE) of image x computed in region R with respect to the reference image r:and the mean of image x for region R,as building blocks.

RMSE(x;r;R)=1|R|∑i∈R(xi−ri)2,

MEAN(x;R)=1|R|∑i∈Rxi,

We required three conditions for an image θ to satisfy the convergence criteria. First, the RMSE of the whole object satisfiesSecond, the background RMSE satisfiesThird, in each VOI Ri the AEM (absolute error of the mean) satisfiesHere r denotes the converged BSREM reference image, W the marginally eroded whole object VOI, B the background VOI and Ri specific VOIs (“tumors,” “spheres,” “white/gray matter,” etc.). Note that metrics therefore only took into account regions of interest in the phantom and excluded the non-radioactive background. The thresholds were chosen to be well-below clinically accepted variation while avoiding being too stringent and depend on numerical errors. In particular, “converged” images were visually very similar, except possibly in the non-radioactive background.

RMSE(θ;r;W)MEAN(r;B)<0.01.

RMSE(θ;r;B)MEAN(r;B)<0.01.

3.6 Infrastructure

To run a public competition, several components are needed: registration forms, template code, data hosting, submission method, standardised and secure GPU-enabled continuous integration (CI), and live leader boards.

We used a dedicated GitHub repository for providing a starter code, and to create private repository forks per participating team. For team registration, GitHub issue forms were used. Hosting everything under one GitHub organisation also allowed secure access to our self-hosted Continuous Integration (CI) runner, which has the required docker image [34, 35], training data, and a GPU. Code uploaded by teams is added to a sequential queue to run on this CI.

All of the reconstructions were run on our standard platform: 32 core CPU, 100 Gb RAM, NVIDIA A100 40 Gb PCIe. The CI machine also runs a web server from which teams could see live TensorBoard logs of their CI runs, and download training data locally. After the end of the submission period, test datasets were added to the machine and all submissions were rerun to produce a final leader board.

All code for this backend machine is available in the corresponding backend repository [36].

4 Results

4.1 Submitted algorithms

Four teams submitted a total of nine algorithms for evaluation, summarised in Table 3. Their submissions are linked in the table which should be referred to for all details. All submissions are gradient-based iterative algorithms but differ in a range of aspects. For the core algorithm, teams used Gradient Descent (GD) with stochastic approximation (Stochastic Gradient Descent; SGD), variance reduction (SAGA, SVRG), as first adapted to PET in [37], and acceleration (Katyusha). See [38] for an overview on these concepts. All teams experimented also with preconditioning ranging from preconditioners motivated from MLEM, diagonal approximations of the RDP Hessian and various combinations of the two. Regarding step-sizes, the teams used a-priori defined decaying step-sizes as are common in subset-based algorithms like BSREM but also experimented with others like Barzilai–Borwein, Armijo line-search and the Distance over Weighted Gradients (DoWG) method [39]. Various teams considered advanced warm-starting options computed by other iterative algorithms as well as deep learning (DL) methods.

Table 3. Overview of submissions sorted by team name in alphabetical order.

The teams also fine-tuned subset-selection methods like the number of subsets and in which order they should be picked. All teams agreed that the number of subsets should be scanner specific and depend on various factors like trading off computational resource and prime factors of the number of views. Typical chosen values are in the range of 20–50. For the actual selection, teams used randomly selected sequences with and without replacement, established deterministic sampling (e.g., Hermann–Meyer) as well as novel variants. It should be noted that some teams also decided to seed the randomness making their stochastic gradient descent algorithm in practice deterministic.

More details on the MaGeZ submission can be found in [40]. The Tomo-Unimib was refined after the challenge ended in [41].

4.2 Algorithms performance and insights

All submitted algorithms, see Table 3, were run until their solution obtained a specific relative error on all selected metrics or their runtime exceeded 1 h, see Section 3.5. Details on ranking method are given in Section 3.4. In Figure 2 we show the rankings of the algorithms with an error bar indicating variability in performance among runs.

The wide variety of the data provided led to varied algorithm performance across different phantom regions. Some algorithms performed well on some datasets, while struggling with others due to poor parameterisation, as can be seen in Figures 3, 4. Notably, while team rankings remained consistent across all runs, there was significant variability in ranking positions among individual submissions from the same team. This highlights the difficulty to maintain performance across such heterogeneous datasets, with potential significative gains to be obtained by accurate parametrisation for each scanner. In general the algorithms from the winning team outperformed all the other algorithms. A selection of metrics for all submissions on 6 different dataset for the whole object and volume of interest is shown in Figure 5, together with the applied threshold.

As discussed in Section 4.1 all submitted algorithms are variants of gradient descent with stochastic approximation, which have acquired popularity recently in literature, [38]. For example, MaGeZ, UCL-EWS and SOS teams all used the SVRG algorithm, nevertheless achieving different performance due to choices in pre-conditioning, subset selections and step sizes.

4.3 Outcomes of the challenge

The outcomes of the challenge were multifaceted. Beyond the competition itself, the event served as a valuable learning experience. One of the key takeaways was the availability of open-source implementations. All three top-scoring teams made their code publicly accessible. While some of these implementations may lack thorough documentation, they are available for exploration and reuse. This openness is a promising step toward future research and development, and some of their contributions are part of the same research topic of Frontiers in Nuclear Medicine. For those looking to get started in this area, these implementations offer a solid foundation. To support reproducibility and further experimentation, we have provided Docker images [34] containing all necessary software.

Another key outcome is the curated phantom data collection which is available on our website, covering a variety of scanners, phantoms and tracers. Additionally, the framework that we developed for the challenge can serve as a starting point for others interested in organising similar challenges.

During the organisation of the challenge, we also encountered and addressed several issues in our own software. This led to significant improvements, including the development of numerous scripts for data ingestion, validation, and visualisation, which may prove useful in many other instances. Major weaknesses in the performance of some algebraic operations, as well as exposure of the underlying data arrays to Python, were identified and resolved after the challenge.

The PETRIC challenge was concluded with a corresponding focused workshop as part of the 2024 IEEE Nuclear Science Symposium and Medical Imaging Conference in Tampa, Florida, USA.¹

5 Discussion and outlook

The PETRIC challenge was successful and achieved several of its aims. The relatively low turnout of participants may be attributed to factors such as the tight timeline from the announcement of the challenge (June 10th 2024) to its start and conclusion in the same year, June 28th and September 30th respectively; but also to the partial availability of training data at the start of the challenge.

One of the aims of this challenge was to evaluate the performance of the reconstruction algorithms on real data. The datasets in this challenge were acquired from several different scanner models and vendors, and different phantoms to cover a wide range of real-world scenarios. The data varied strongly in terms of imaging parameters (e.g., time-of-flight vs. non-time-of-flight) and noise (e.g., due to different scan durations and phantom fillings). This posed a significant challenge for all the submitted algorithms because they had to perform well for a very heterogeneous set of data.

This allowed us to assess how well the performance of the different algorithms generalises to different applications. As mentioned above, to ensure the challenge was as fair as possible, data used for the final scoring was obtained after the submission deadline of the challenge. This was only possible due to the generous contribution of data from several different research sites.

One major difficulty with using data from different scanners was to convert the vendor data to a common open data format which is compatible with SIRF and STIR, respectively. In this occasion we used STIR interfile format for PET projection data. Some vendors and scanner models required closed-source tools only available based on research agreements with the vendors and detailed knowledge about the vendor data formats. This restricted the number of research sites which could contribute datasets to the challenge. In the future more scanners may become available thanks to the new initiative on emission tomography standardisation, of which one aim is to design a new file format, the PET Standardisation Initiative’s Raw Data (PETSIRD) [42].

Using real data meant that no ground truth reference solution was available and therefore the scoring of algorithms needed careful design. Rather than setting up a challenge which looks for the algorithm yielding the best possible image, we designed the challenge to find the algorithm which reaches a known solution in the shortest possible computational time (within an existing reconstruction framework). Reconstruction speed is also an important quality metric of how useful an image reconstruction algorithm is for clinical applications. Of course, in addition to algorithm design, software implementation of projectors and numerical operations have large impact on reconstruction speed. This would be an interesting aspect for future challenges, but would need careful design of the evaluation criteria as well as software infrastructure.

Use of numerical simulations instead of real data would have provided ground truth solutions. Nevertheless, we focused on the assessment of the performance of the algorithms on real data which makes the translational step from scientific development to application on clinically relevant data smaller.

The focus of this challenge was on classical optimisation schemes to solve a well-defined optimisation problem. This allowed us to evaluate the different submissions even without having a ground truth solution. Nevertheless, this made this challenge less suited for deep-learning-based image reconstruction algorithms. While implementing learned reconstruction methods in the STIR, SIRF and CIL software is already feasible [43], current software limitations affect run-time. An important aim of our community is to extend STIR, SIRF and CIL to support memory-efficient integration with Python CUDA libraries such as PyTorch, before organising a challenge which is also suitable for deep-learning-based submissions. Such a challenge would likely concentrate on metrics evaluating image quality (as opposed to convergence), including lesion detectability, contrast and noise. Due to the nature of deep-learning methods, data from most “standard” phantoms is not suitable, and will likely have to be replaced by a combination of simulated and measured clinical data, with associated difficulties in data access.

Given the successes of PETRIC, we have launched PETRIC2, which is open at the time of writing and is similar to PETRIC but uses data at reduced count level. Additional challenges are also considered, for example by adapting the methodology for SPECT reconstruction.

6 Conclusion

In this article, we presented PETRIC, the first challenge for PET image reconstruction using measured PET datasets from a range of PET scanners. Its solid foundations allow further adaptation for different types of evaluation of the performance of image reconstruction methods on existing or new datasets from any PET scanners, radionuclides, and phantoms.

Generative AI statement

The author(s) declared that Generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

Edited by: Georgios Angelis, The University of Sydney, Australia

Reviewed by: Alan Miranda, University of Antwerp, Belgium

Evangelos Papoutsellis, The University of Manchester, United Kingdom

Footnotes

¹Recordings of the presentations at the workshop, including of the three top-ranked teams, are available at https://www.ccpsynerbi.ac.uk/petric-workshop-2024-recordings/.