Concept development of an on-chip PET system

Monte Carlo simulation

We model the interaction of the proposed system with a Fluorine-18 positron source in an MCS built with the Geant4 Application for Emission Tomography (GATE) tool [16,17,18]. GATE enables the creation of MCSs in the field of nuclear medicine through a macro language that controls the experimental settings. It is built as a wrapper around the Geant4 simulation toolkit that enables simulating ”the passage of particles through matter” for a wide range of physics processes, particles, and materials over a broad energy spectrum [19,20,21].

A GATE simulation consists of the following parts that are described in more detail in the next sections: a scanner geometry, a phantom, material properties, physics processes, boundary surfaces, and a primary particle source.

Geometry

In GATE, the concept of a system plays a crucial role if one wants to store information about particles and physical processes in the simulation. A system is defined as a template for predefined scanner types with specific geometries. Different geometrical shapes are organized in a tree-level structure for these scanners to build up the final geometry.

In this work, we use the most generic system in GATE, the scanner. In our case, the scanner is defined as a box-shaped volume and placed in the world volume. The purpose of the scanner volume is to encapsulate our proposed PET system. The box-shaped detector volume is placed inside the scanner volume and repeated four times with a ring repeater around the z-axis. Inside the detector volume, the box-shaped crystal volume is placed and repeated two times with a linear repeater. Table 1 contains the lengths, translations, and materials of each volume. We evaluate two different crystal thicknesses, 13 mm and 26 mm. All dimensions of the crystals are chosen such that arrays of commercially available SiPMs fit on the surfaces. The properties of the respective materials of the different volumes are described in the next section.

Table 1 Geometry setup of the GATE simulation

In GATE, it is important to attach the created volumes to the system to be able to record hits in them. In our case, the detector volume is attached to the first level of the scanner system and the crystal volume to the second one. To both of these volumes, we attach a crystal sensitive detector that records the hits in them.

Figure 1 depicts the simulation setup viewed from the front and the side.

Fig. 1figure 1

Simulation setup viewed from the front and the side. The yellow volumes represent the crystals, the green ones the epoxy layer around the crystal, and the blue one the volume from which the source position is sampled to create the training dataset

Phantoms

We use different types of phantoms depending on the purpose of the dataset that is created by the simulation.

For training, we use a simplified model of an OOC device consisting of a box-shaped water volume surrounded by a PMMA box that is placed in the middle of the detector. The PMMA box has side lengths of 10 mm \(\times\) 26 mmm \(\times\) 76 mm while those of the water box are 4 mm shorter each.

For evaluating the sensitivity of the system, we instead add a 1 mm thick and 104 mm long water cylinder to the simulation that is placed along the z-axis of the scanner.

To demonstrate that our system is able to capture the small volumes found in OOC devices, we create a more complex OOC phantom consisting of the PMMA box mentioned above that contains four water spheres with a radius of 2 mm each, which are distributed along the z-axis.

In “Sources” section, the sources that are combined with each of the phantoms are described.

Materials

The materials of the world and scanner volumes are set to vacuum, the detector’s to epoxy, and the crystal’s to LYSO. In this way, one detector consists of two LYSO crystals that are surrounded by a 0.1 mm thick layer of epoxy, which acts as the glue and optical medium between the crystals and SiPMs.

In Table 2, the material properties are described, and in Table 3, the scintillation properties of LYSO are shown.

Table 2 Properties of the materials used in the GATE simulationTable 3 Scintillation properties of LYSOPhysics and cuts

As a physics list, the Electromagnetics (EM) constructor with option four is chosen, which uses the most accurate standard and low-energy models available in Geant4. Table 4 shows the enabled physics processes with their selected models. We set a cut of 0.1 mm for gammas, electrons, and positrons in the scanner volume.

Table 4 Added processes in the GATE simulationSurfaces

We use Geant4’s unified model to define the surfaces in the simulation. The surfaces between the detector and crystal volumes are dielectric-dielectric ones with a ground finish and a sigmaalpha value of 0.01 corresponding to a typical polished crystal. Their specular lobe constant is set to 1.0.

To detect optical photons in GATE, it is necessary to use dielectric-metal boundaries. As we want to detect the optical photons that are leaving the Epoxy layer and thus entering the SiPMs, we add dielectric-metal surfaces to the boundaries between the scanner and detector volumes. These surfaces have specular lobe, specular spike, and backscatter constants of 0.0, a reflectivity of 0.0, and an efficiency of 1.0.

Sources

To create the training and testing datasets, we add a source that emits Fluorine-18 positrons with an activity of 1000 Bq to the simulation setup. The position of the source is sampled from the box-shaped water volume described in “Phantoms” section.

To evaluate the spatial resolution of the system, we create another dataset where 21 F18-point sources with a total activity of 1000 Bq are arranged in a 7 \(\times\) 3 grid with a distance of 10 mm in between each source, which are placed in the box-shaped water volume.

A 104 mm long line source is placed along the z-axis inside the water cylinder to determine the sensitivity of the system.

To create the more complex OOC phantom, we place four hot sphere-shaped sources with radii of 0.4 mm, 0.5 mm, 0.6 mm, and 0.7 mm inside the four water spheres described in “Phantoms” section. These hot sources each have an activity concentration of 1000 Bq/mm\(^3\), and are surrounded by cold sphere-shaped sources with a radius of 2 mm and an activity concentration of 100 Bq/mm\(^3\).

Dataset creation

In the first step, we run the GATE simulation described in the previous section in parallel. In the second step, we post-process the hits output files from the GATE simulation runs. For each primary event, in which one positron is emitted, zero to two samples of the training dataset are created. The amount of created samples depends on the number of detectors in which the gamma-rays interact. If there is a scintillation event in at least one detector, the position of the first interaction of the gamma-ray as well as the corresponding light patterns that emerge on the surfaces of the detectors, are saved. To create the light patterns, we simulate the SiPMs that are attached to the surfaces of the detector. This procedure is described in more detail in the next section. Table 5 gives an overview of the created datasets. Depending on how many different datasets are needed to evaluate a design choice, we use either the smaller (100k) or larger (1M) versions.

Table 5 Overview of created datasetsLight pattern creation

To create the light patterns that emerge on the surfaces of the detector, we simulate the SiPMs that cover the surfaces of the detector with SimSiPM, a C++ library for SiPM simulation [22]. We want to evaluate the performance of our proposed system with three different options for the size of the photosensitive area of the SiPMs. The SiPMs simulated in this work come from Hamamatsu’s S14161 series with photosensitive area sizes of either 3, 4, or 6 mm [23]. From now on, we refer to the SiPMs as small (3 mm), medium (4 mm), and large (6 mm). The properties of the SiPMs are described in Table 6.

Table 6 Properties of the simulated Hamamatsu SiPMs [23]

The light patterns are created with the following procedure:

1.

Assign the corresponding SiPM to each photon leaving the epoxy layer depending on the photon’s position.

2.

Add photon times for each SiPM on each surface.

3.

Integrate over the SiPM signal to compute the output value of each SiPM on each surface.

4.

Pad the resulting light patterns with zeros such that they all have the same shape. This step is necessary to be able to use the stacked light patterns as input to the CNN.

Six light patterns recorded with different SiPMs sizes of 13 mm and 26 mm thick crystals are depicted in Fig. 2.

Fig. 2figure 2

Light patterns recorded with different SiPMs sizes of 13 mm thick crystals (left) and 26 mm thick crystals (right). The light patterns recorded with the SiPMs on all surfaces of the detector except the front one are shown. The top light patterns are created with a photosensitive area size of the SiPMs of 3 mm, the middle ones with 4 mm, and the bottom ones with 6 mm. The light patterns are padded such that they all have the same square shape. The red dot is the scintillation position of the gamma-ray inside the crystal projected onto each surface

Scintillation position prediction

With the created dataset described in the previous section, we train a CNN that predicts the gamma-ray interaction position inside the crystal. The input to the network is the stacked light patterns recorded with SiPMs on the detectors’ surfaces. With the light pattern images as input, the network should predict the gamma-ray interaction positions.

Baseline method

In addition to our deep learning-based approach, we implement a simple centroiding-based method to serve as a baseline.

For every sample in the dataset, the baseline method performs the following steps to determine the scintillation position:

1.

Compute the centroids of every light pattern with image moments.

2.

Take the mean of the centroids of corresponding light patterns (top-bottom, left-right) to compute x-z and x-y centroids. This step is not applicable for the back light pattern as there is no corresponding front light pattern.

3.

Take the mean of the corresponding dimensions of the x-y, x-z, and y-z centroids to compute the x-, y-, and z-centroid.

4.

Stack the x-, y-, and z-centroid to create the predicted x-y-z position.

The predicted position is then compared to the ground truth value for every sample in the dataset to compute the MAE value of the baseline method.

Training pipeline

We implement the CNN training pipeline with PyTorch Lightning [24]. The code is organized with subclasses of the LightningDataModule and the LightningModule. In the LightningDataModule subclass, the train and test datasets are loaded, the train dataset is split into a train and a validation dataset, and the PyTorch DataLoaders are set up. In the LightningModule subclass, the network architecture, loss function, and optimizer are initialized. Additionally, the training, validation, and test steps are defined, including logging the loss and metrics. Weights & Biases [25] is used to create experiment sweeps and log results. The training pipeline of the scintillation position prediction is depicted in Fig. 3.

Fig. 3figure 3

Training pipeline of the scintillation position prediction. The light patterns recorded on each surface are stacked along the channels dimension and used as input to the CNN. The CNN should predict the gamma-ray interaction position from the stacked light patterns. The MAE loss is computed with the predicted and ground truth positions and backpropagated through the network to update its weights

Experiment setup

In this work, the design of our proposed On-Chip PET system is optimized. We investigate the influences of the network architecture, crystal thickness, SiPM size, and the number of surfaces covered with SiPMs on the scintillation-position prediction performance. For the network architecture, we selected variants from the ConvNeXt [26], ResNet [27], and EfficientNet [28] families with different depths, as these architectures have shown great performances on image data. The parameters given in Table 7 are used for every experiment run. The experiments are performed on a machine with four NVIDIA GeForce RTX 3090 GPUs where one training takes approximately one hour to run.

Table 7 Parameters used for every experiment runSensitivity

We use the line source along the z-axis described in “Sources” section to determine the sensitivity of the system. The source has an activity of 1000 Bq and is simulated for 10 s. The sensitivity is computed by dividing the number of coincidence events where two interactions are recorded in two different detectors by the total number of events.

Reconstruction

We evaluate the reconstruction performance of our proposed system with a grid of point sources and the more complex OOC phantom, as described in “Phantoms” and “Sources” sections.

For the reconstruction, it is necessary to be able to create LORs, which is only possible if the back-to-back gamma-rays interact in two detectors at the same time. For each of those sample pairs, the two corresponding scintillation positions are predicted with a trained network. All pairs of predicted scintillation positions are stored for further processing.

We perform the following steps for SART, an algebraic reconstruction method that shows good reconstruction performance in cases where there is limited data available [14]:

1.

Load all predicted pairs of scintillation positions from disk.

2.

Compute the distance and angle of the LOR to the origin for each LOR that is defined by the pair of predicted scintillation positions.

3.

Create the sinogram from the LORs by computing the 2D histogram of the distances and angles with a bin size of 400 in both dimensions.

4.

Generate the corresponding reconstructed image by running five iterations of SART implemented in scikit-image [30].

留言 (0)

沒有登入
gif