PIPET: A Pipeline to Generate PET Phantom Datasets for Reconstruction Based on Convolutional Neural Network Training

Abstract

There has been a strong interest in using neural networks to solve several tasks in PET medical imaging. One of the main problems faced when using neural networks is the quality, quantity, and availability of data to train the algorithms. In order to address this issue, we have developed a pipeline that enables the generation of voxelized synthetic PET phantoms, simulates the acquisition of a PET scan, and reconstructs the image from the simulated data. In order to achieve these results, several pieces of software are used in the different steps of the pipeline. This pipeline solves the problem of generating diverse PET datasets and images of high quality for different types of phantoms and configurations. The data obtained from this pipeline can be used to train convolutional neural networks for PET reconstruction.

1. Introduction

In nuclear medicine, there is a functional medical imaging technique called Positron Emission Tomography (PET). This functional medical imaging technique is particularly useful in obtaining information from biological (or functional) processes, such as metabolism, blood flow, absorption of pharmaceuticals, etc.

In order to obtain biological information, we need a compound that will be injected into the patient’s body. These compounds are referred to as radiotracers, chemical compounds in which one or several atoms are replaced by a radioisotope, an element with an unstable nucleus. In the case of PET, the most widely used radiotracer is fluorodeoxyglucose F18 (FDG) [1,2].

This molecule is analogous to a glucose molecule with the replacement of a hydroxyl radical with the radionuclide fluorine-18, which emits positrons through a ${\beta }^{+}$ decay. This type of functional imaging is extensively used in oncology applications, as tumors or cancerous cells consume vast quantities of glucose to enable their uncontrolled reproduction. Therefore, when a patient receives a radiotracer dose, it is usually detected in the area where the cancerous cells are located. As a consequence of this phenomenon, those areas are easily detected by a PET scanner as regions with high activity concentrations, also known as hot regions, due to the radioactivity of the radiotracer.

The physical process that PET relies on is the detection of two 511 KeV photons generated from the matter–antimatter annihilation of an electron with a positron. In the data acquisition step of a PET scanner, those photons are detected at an angle at 180° forming a line called the Line Of Response (LOR), which defines the path followed by the pair of detected photons emitted in the annihilation event.

In order to calibrate a PET scanner, a scan with an imaging phantom (or just phantom) is needed. An imaging phantom is a device specifically designed to simulate a standardized human body or part of it and imitate the response of human tissues and organs when injected with FDG. Therefore, imaging phantoms are used to analyze, evaluate, and tune PET scanners. To that effect, a series of tests, of which large part have to be performed on the reconstructed image, are specified by the NEMA [3,4,5] protocols.

The raw data obtained from a PET scan must be processed to obtain the required reconstructed image. Different methods exist to perform this reconstruction. There are simple, non-iterative methods, such as Inverse Radon Transform [6] or Filtered Back Projection (FBP) [7]. These methods perform reasonably well, but the results are usually not of high quality. There are also more complex iterative reconstruction methods, such as Maximum Likelihood Expectation Maximization (MLEM) [8] or Ordered Subset Expectation Maximization (OSEM) [9]. These methods are usually computationally expensive but result in high-quality reconstructions of the patient or phantom used in scan acquisition.

The rapid and successful advancements in the field of Artificial Intelligence (AI) that have been enabled by the use of neural networks (NNs) have also been adopted in the medical imaging field. In state-of-the-art PET imaging reconstruction, the use of convolutional neural networks (CNNs) has delivered good results in the reconstruction of images from data acquired with PET scans as input for the neural network [10,11].

The use of neural networks has some caveats: they are prone to overfitting, they are black-box models, and they require a large number of data during the training step. In the case of PET image reconstruction, a large quantity of images is needed in order to train CNNs to perform the reconstruction task.

There are two ways to acquire PET data and images for studies. The first involves the acquisition of real scans from clinical cases, PET data repositories, or data from experimental scanners. In the case of clinical cases or experimental scanners, this requires collaboration with medical centers or research groups and has an associated cost in terms of equipment and work hours of medical professionals. Furthermore, this provides a limited number of data, which limits the studies that could be performed. In the case of PET imaging, while there are multiple repositories with reconstructed PET images, such as The Cancer Imaging Archive [12] or the OpenNeuro database, we have not found any available open repositories with raw data from real PET scans, as commercial scanners usually output the final reconstructed image.

The second method is to simulate PET data acquisition by using different pieces of software from the generation of the phantom to the acquisition and the reconstruction. This method has been shown to provide results which are usually of sufficient quality and comparable to the acquisition of real measurements [13,14,15,16,17]. This method has the disadvantage of requiring vast computational resources to obtain the data. However, it can provide a large number of data with enough computation time.

To address the problem of limited datasets, we intend to create a pipeline that will programmatically generate different types of medical phantom raw data and images with enough diversity to generate a dataset that can be used in future works to train a neural network reconstructor of PET images. This paper is an extended version of our conference paper [18].

2. Materials and Methods

2.1. Pipeline and Phantom Definition

Due to the nature of the different pieces of software that are employed to generate PET phantoms, we will first identify the inputs and outputs of each of the different programs that are used to achieve this goal.

The Figure 1 is a diagram of the pipeline developed.

The steps followed in the pipeline are described below:

• Select the type and the number of phantoms to generate.
• Generate the 3D voxelized activity map and attenuation map for the phantoms; the tuning of the parameters for each phantom are selected at random within the ranges defined for each phantom during generation.
• Move the generated phantoms to a CPU cluster to run the simulations for each generated phantom.
• For each phantom, join the data files of the simulation, extract the global coordinates, correct the DOI, and transform the data into list mode.
• Convert the list-mode data into CASToR format and reconstruct them by using iterative methods.

In PET imaging, the most widely used phantoms are those reported in

Table 1. Types of phantoms and geometries.

Phantom Name	Main Geometry	Secondary Geometry	Number of Secondary Geometries
NEMA	Cylinder	Cylinder	6
Jaszczak	Cylinder	Sphere	6
Derenzo	Cylinder	Cylinder	116
Shepp–Logan	Ellipse	Ellipse	10

In the current implementation, the pipeline must be run individually for each phantom type. Therefore, the code must be run independently to generate the desired number of NEMA, Jasczak, Derenzo, and Shepp–Logan phantoms.

.

Some samples of the mentioned phantoms can be seen in Figure 2.

As shown in the figure above, 3D phantoms have different shapes, sizes, and numbers of secondary geometries, such as cylinders, spheres, or ellipses. Now that we have a clearer view of the different elements that are part of the pipeline developed, we will address each part in detail.

2.2. Phantom Generation

The first step in PIPET involves the generation of voxelized 3D phantoms. In order to achieve this goal, a series of Python [19] scripts were developed to allow for the tuning of the different parameters of these phantoms. In order to generate a diverse dataset for later use as training data in a neural network, it is important to have enough variation in the phantoms to prevent the network from memorizing the patterns, which could lead to overfitting. It is necessary to have similar data but with enough differences that the network tends to avoid learning the data and suffering from overfitting [20].

Considering this fact, it was decided to generate a dataset containing phantoms similar to the ones used in real scan acquisition. As standard medical phantoms’ design and dimensions are well established and defined, the generated phantoms will be referred to as like phantoms in reference to the baseline phantom used. In order to generate the voxelized phantoms, we modified the main geometry size, shape, and positions of the secondary geometries within each generated phantom. The following subsections will explain the particularities and processes implemented to generate them.

2.2.1. NEMA-like Phantom

The NEMA NU 4 [3] is a medical imaging phantom used to calibrate and compare the performance of small-animal PET scanners. This phantom’s main geometry is a cylinder; the secondary geometries are six smaller cylinders of various sizes arranged in a radial disposition from the center of the phantom. The parameters that can be evaluated with this particular phantom are spatial resolution, scatter fraction, count rate performance, and sensitivity, among others.

The parameters shown below are considered to generate NEMA NU 4-like phantoms.

In the standard NEMA phantom, the secondary cylinders inserted into the main cylinder must be equidistant. As a consequence of this, in order to develop a NEMA NU 4-like phantom in which the number of secondary cylinders can be modified, the main cylinder needs to be split into sectors, and the maximum radius of the secondary geometry has to be calculated in order to fit into these sectors. The procedure to generate a NEMA-like phantom is the following:

• Randomly select the main cylinder’s parameters and the number of secondary cylinders based on the ranges defined in

  Table 2. NEMA parameters.

  Phantom Name	Main Geometry	Parameter	Tunable	Range	Secondary Geometry	Parameter	Tunable	Range
  NEMA	Cylinder	Radius (mm)	Yes	30–60	Cylinder	Radius (mm)	Yes	Depends on number of geometries
  		Height (mm)	Yes	20–50		Height (mm)	No
  						Number	Yes	3–8

  .
• Split the main cylinder into sectors, considering the number of secondary cylinders. In addition, calculate the maximum radius contained in a circular sector.
• Randomly select the parameters of the secondary cylinders based on the ranges defined in Table 2.
• Create an empty array with the dimensions of the generated phantom image.
• Iterate through the array considering the voxelized coordinates of the main and secondary cylinders to assign an activity value, assigning the highest activity values to the secondary cylinders.

2.2.2. Jaszczak-like Phantom

The Jaszczak phantom is a medical imaging phantom used to evaluate spatial resolution, image uniformity, and attenuation correction accuracy. The main geometry of this phantom is a cylinder, and the secondary geometries are six spheres of various radii, arranged in a radial disposition similar to the NEMA NU 4 phantom.

The parameters shown below are tuned during the generation of Jaszczak-like phantoms.

The procedure to generate a Jaszczak-like phantom is the following:

• Randomly select the main cylinder’s parameters and the number of secondary spheres based on the ranges defined in

  Table 3. Jaszczak parameters.

  Phantom Name	Main Geometry	Parameter	Tunable	Range	Secondary Geometry	Parameter	Tunable	Range
  Jaszczak	Cylinder	Radius (mm)	Yes	30–60	Sphere	Radius (mm)	Yes	Depends on number of geometries
  		Height (mm)	Yes	20–50		Height (mm)	No
  						Number	Yes	3–8

  .
• Split the main cylinder into sectors, considering the number of secondary spheres. In addition, calculate the maximum radius contained in a circular sector.
• Randomly select the parameters of the secondary cylinders based on the ranges defined in Table 3.
• Create an empty array with the dimensions of the generated phantom image.
• Iterate through the array considering the voxelized coordinates of the main cylinder and secondary spheres to assign an activity value, assigning the highest activity values to the spheres.

2.2.3. Derenzo-like Phantom

The primary geometry of the Derenzo phantom is a cylinder, with secondary geometry consisting of smaller cylinders arranged in a triangular pattern. In our Derenzo-like phantoms, we use triangular numbers to calculate the number of cylinders to include in each sector according to the following formula:

$$T_n\equiv \sum _{k=1}^n={\displaystyle \frac{1}{2}}n(n+1)=\left({\displaystyle \genfrac{}{}{0pt}{}{n+1}{k}}\right)$$

(1)

where $\left({\displaystyle \genfrac{}{}{0pt}{}{n}{k}}\right)$ is a binomial coefficient.

This phantom is used to evaluate the imaging system resolution by measuring the full width at half maximum (FWHM) between the secondary cylinders in each sector, as within each sector, the series of cylinders are equidistant and have the same radius. The following parameters are tuned during the generation of Derenzo-like phantoms.

The procedure to generate a Derenzo-like phantom is the following:

• Randomly select the main cylinder’s parameters and the number of sectors into which to split the main cylinder based on the ranges defined in

  Table 4. Derenzo parameters.

  Phantom Name	Main Geometry	Parameter	Tunable	Range	Secondary Geometry	Parameter	Tunable	Range
  Derenzo	Cylinder	Radius (mm)	Yes	40–60	Cylinder	Radius (mm)	Yes	Depends on number of geometries
  		Height (mm)	Yes	20–50		Height (mm)	No
  						Number	Yes	3–7

  .
• Calculate the triangular number progression, considering the number of sectors, and assign the value to the number of cylinders that each sector will contain.
• Assign the radius of each cylinder, considering the number of cylinders to fit in each sector.
• Create an empty array with the dimensions of the generated phantom image.
• Iterate through the array considering the voxelized coordinates of the main cylinder and secondary spheres to assign an activity value, assigning the highest activity values to the cylinders.

2.2.4. Shepp–Logan-like Phantom

The Shepp–Logan phantom is a mathematically defined phantom that consists of a series of ellipses, existing only in simulations. This phantom is used to characterize the detector’s spatial resolution, noise characteristics, artifact generation, and contrast resolution. The main geometry is an ellipse, and the secondary geometries are a set of ellipses. It resembles the anatomical structures of the human torso, considering organs such as the lungs, heart, or liver.

The following parameters are tuned during the generation of Shepp–Logan-like phantoms.

The major and minor axis of the phantom are selected at random. Therefore, these ranges are not been included in the above table. The main condition is that the ellipses do not exceed the dimensions of the voxelized image.

The procedure to generate a Shepp–Logan-like phantom is the following:

• Define the major and minor axis ranges for each ellipse, avoiding overlap between the different secondary geometry ellipses based on the ranges defined in

  Table 5. Shepp–Logan parameters.

  Phantom Name	Main Geometry	Parameter	Tunable	Secondary Geometry	Parameter	Tunable
  Shepp–Logan	Ellipse	Major axis	Yes	Ellipse	Major axis	Yes
  		Minor axis	Yes		Minor axis	Yes
  					Number	Yes

  .
• Calculate the voxelized coordinates of each ellipse.
• Create an empty array with the dimensions of the generated phantom image.
• Iterate through the array considering the voxelized coordinates of the main ellipse and secondary ellipses to assign an activity value.

2.3. PET Acquisition Simulation

In this pipeline step, the generated like phantoms were used to simulate PET acquisition by using the GEANT4 Application for Tomographic Emission (GATE) V-9.2 [21] software suite, which relies on numerical Monte Carlo methods. In this instance, data acquisition was run for a total time of 1500 s for each like phantom, using the scanner described in the article [22]. This scanner was designed to obtain breast scans for cancer screening; the FOV was 200 mm in the transaxial dimensions and 50 mm in the axial one, and the phantoms from the previous step were generated with this detector in mind.

In Figure 3, we can see the simulated scanner in GATE with a generic cylindrical phantom.

The data obtained from these simulations were stored in a ROOT file that registered the coordinates, energies, and many other variables of PET acquisition. ROOT [23] is a framework developed by CERN (Conseil Européen pour la Recherche Nucléaire) for scientific data analysis, used mainly in the nuclear physics field.

2.3.1. Data Processing

Once the data from PET acquisition were obtained, they were processed by using the ROOT framework to extract the X, Y, and Z coordinates, energies, and time for each coincidence event in the simulation. In addition, the Depth of Interaction (DOI) of these data was corrected. The DOI is an effect of the detection of the pair of photons that conforms to each coincidence event, as the depth (Z coordinate) in the scintillator crystal affects the coordinates for the oblique coincidence event detected. This correction is necessary to improve the spatial resolution of the images obtained. After processing, the data were stored in list mode; this data storage method orders and aggregates the information of each pair of photons detected in each coincidence event and is a requirement for the reconstruction step of the pipeline.

2.3.2. Reconstruction

In the pipeline’s last step, the list-mode data had to be transformed into a format recognized by CASToR. The list-mode data were reconstructed by CASToR [24], tomographic reconstruction software, by using the MLEM iterative method. Twenty iterations were used, and both attenuation and sensitivity corrections were applied. The final images were of dimensions 200 × 200 × 50 mm³, with a resolution of 1 voxel/mm³

In order to be able to use these reconstructed images to train a 3D or 2D reconstruction convolutional neural network, further processing would have to be applied. However, this task is out of the scope of this work.

2.3.3. Evaluation

Several methodologies were employed to evaluate the voxelized images against the reconstructions. A visual check of the resemblance between the voxelized images and the reconstruction was performed to compare the activity distribution between the two cases. Additionally, a line profile was also obtained. A line profile shows the intensity of an image along a specified line or path. We compared the line profile of both the voxelized and reconstructed images to validate the results from the visual check.

In order to obtain a quantitative metric for all available images, another method of evaluation was the use of the Structural Similarity Index (SSIM), which indicates the similitude between two images. Its range is between 0 and 1, indicating that the maximum value equals a perfect copy of the reference.

$$SSIM(x,y)={\displaystyle \frac{(2{\mu }_x{\mu }_y+c_1)(2{\sigma }_{xy}+c2)}{({\mu }_x^2+{\mu }_y^2+c_1)({\sigma }_x^2+{\sigma }_y^2+c_2)}}$$

(2)

where the following apply:

• ${\mu }_x$ is the pixel sample mean of x.
• ${\mu }_y$ is the pixel sample mean of y.
• ${\sigma }_x^2$ is the variance of x.
• ${\sigma }_y^2$ is the variance of y.
• ${\sigma }_{xy}$ is the covariance of x and y.
• $c_1={\left(k_{1}L\right)}^2$ and $c_2={\left(k_{2}L\right)}^2$ are two variables to stabilize the division with a weak denominator.
• L is the dynamic range of the pixel values (typically, this is ${\displaystyle 2^{\#bitsperpixel}-1}\left)\right)$.
• $k_1=0.01$, and $k_2=0.03$ by default.

Lastly, we also calculated the Fréchet Inception Distance (FID); this metric is used to check the quality and diversity of datasets, as it performs feature extraction and utilizes the means and covariance matrix of the feature matrix of each image. For the interpretation of this metric, the closer it is to zero, the more similar the images from the datasets are. When the FID < 10, the datasets are highly similar and of high quality [25,26].

$$FID={|\mu -{\mu }_w|}^2+tr(\Sigma +{\Sigma }_w-2{\left(\Sigma {\Sigma }_w\right)}^{{\displaystyle \frac{1}{2}}})$$

3. Results

Once PIPET was implemented, it was decided to simulate a total of 400 phantoms. These included 100 NEMA-like, 100 Jasczack-like, 100 Derenzo-like, and 100 Shepp–Logan-like phantoms, as these different phantoms can provide enough variability and a diverse dataset to be used in future works as the training dataset to implement a CNN reconstructor. Given the large number of scan acquisition processes to simulate, it was decided to run the GATE simulations on a CPU cluster. Since these processes are single-threaded, performance scales with the number of available CPU nodes and threads. The more computational resources can be allocated, the more efficient the simulations are, as each thread handles a separate task. This approach ensures that the workload is distributed across multiple CPUs, significantly reducing the overall simulation time. By leveraging the power of parallel processing, the CPU cluster maximizes efficiency, making it ideal for handling the computational demands of GATE simulations. The total computation time required to obtain the data was equivalent to approximately 9000 h, with the most significant computational expense of the simulation being PET acquisition with Monte Carlo methods [27,28], which took nearly 8600 h; then, data processing lasted approximately 100 h, and the reconstructions took approximately 270 h. It is well known that Monte Carlo methods are slow and require a large quantity of resources for computation.

In the Figure 4, we can see some of the results obtained in the first step of the pipeline, the generation of the voxelized 3D phantom, and last step, the reconstruction in CaSTOR.

Additionally, line profiles were obtained for the samples shown in Figure 5. The line profiles were obtained for each phantom in the region of higher intensity.

The SSIM was computed for all the phantoms; then, the mean for each type of phantom was calculated. The results of the SSIM metric are reported in

Table 6. Mean Structural Similarity Index for each type of phantom.

Phantom	Mean SSIM
NEMA	0.851
Jaczaczk	0.934
Derenzo	0.879
Shepp–Logan	0.487

.

The FID was computed for each type of phantom; as this metric must be used with datasets, we used each type of phantom as a different dataset. The results obtained are reported in

Table 7. Fréchet Inception Distance for each type of phantom.

Phantom	FID Score
NEMA	0.056
Jaczaczk	0.484
Derenzo	0.061
Shepp–Logan	0.529

.

4. Discussion

By comparing the 3D reconstructions of the simulated phantoms with the activity maps used in the simulation, as shown in Figure 4, we can see how both show the same distribution of activity. In order to validate the results obtained, a manual check was performed on all the generated phantoms to confirm that the voxelized images and their reconstructions were similar and that the radiotracer was found in the same regions for all the cases. Qualitatively, it was observed that for all the images, the voxelized phantoms and the reconstructed images were identical in all cases. Due to the nature of the PET acquisition simulation, the reconstructed images were of lower quality than the voxelized images.

Figure 5 shows that the NEMA, Jaczaczk, and Derenzo reconstructed images showed good agreement with the original voxelized images because the intensities in the regions of interest were similar. In the case of the Shepp–Logan line profile, the line profiles were quite different; this was due to the more complex geometry of the Shepp–Logan phantom; the acquisition was more computationally expensive, which means we had lower total coincidence for the same simulation time compared with the NEMA, Jaczaczk, and Derenzo cases, resulting in a worse line profile. Nonetheless, the region of maximum intensity showed fairly good agreement in the voxelized and reconstructed phantoms. The acquisition times should be longer for the Shepp–Logan-like phantoms to address this issue. From the information provided from the line profiles, we can confirm that the reconstruction obtained from the simulated data in GATE was successful, as the direct comparison of the voxelized image and the reconstruction was consistent.

In regards to the quantitative metrics employed, from the SSIM, it can be interpreted that the voxelized and reconstructed images were similar in the NEMA, Jaczaczk, and Derenzo cases. The SSIM result for the Shepp–Logan phantom was lower than the other three due to the lower coincidence counts, as we have already discussed. A longer acquisition time is required to obtain similar results. In general, the values are high enough to confirm the agreement between the voxelized and reconstructed images. The values of the SSIM indicate that for the generated phantoms, the voxelized and reconstructed images were in agreement.

Finally, from the FID, it can be interpreted that in the NEMA and Derenzo cases, the voxelized images and the reconstructions were quite similar, as the FID score value was close to zero. For the Jaczaczk and Shepp–Logan cases, the FID score was also relatively low. However, it indicated that the reconstructions of these phantoms were not as good as in the other cases. Nonetheless, considering that the results were below the value where they would be considered similar and of high quality, it can be stated that this metric is also in accordance with the results obtained so far.

In order to thoroughly validate the results obtained with the pipeline, the ideal criterion would be to compare the simulated dataset with data from real PET scans. Unfortunately, due to limitations in the availability of PET scans and the impossibility of obtaining data from clinical sources due to data privacy concerns, the data, at this time, could not be compared against real data. In order to validate the results, several phantoms would need to be built to test them properly. However, according to research conducted by other authors, a reconstruction neural network can be trained for reconstruction tasks by using synthetically generated datasets [29,30].

Therefore, it can be concluded that from the results of the visual inspection, line profiles, SSIM, and FID, the pipeline works accordingly with the intention of its development and it provides a reliable and efficient way of obtaining synthetic PET data with sufficient variability, producing a diverse and large dataset to be used in the training of a neural network designed with the task of reconstructing PET images.

Further processing will be applied to the generated dataset for use in future works, as state-of-the-art CNN PET reconstructors work with 2D [31,32] slices and thus are able to work with the raw data obtained so far; rebinning algorithms will need to be applied to the data, as the addition of events with some degree of obliquity on the transaxial plane will have to be considered.

5. Conclusions

PIPET is a pipeline for creating synthetic 3D PET data. It implements the programmatic generation of voxelized like phantoms, PET acquisition simulation using Monte Carlo techniques, and the reconstruction of PET images using well-known iterative algorithms. PIPET enables the generation of an indefinite number of PET phantoms, combining several pieces of software used in this field. The resulting phantoms are used to generate a dataset that includes PET acquisition raw data and PET reconstructed images that can be used as the training dataset for neural networks with the task of reconstructing PET images.