Optimization of the Algorithm for the Implementation of Point Spread Function in the 3D-OSEM Reconstruction Algorithm Based on the List-Mode Micro PET Data

Abstract

Positron emission tomography (PET) is a popular research topic. People are becoming more interested in PET images as they become more widely available. However, the partial volume effect (PVE) in PET images remains one of the most influential factors causing the resolution of PET images to degrade. It is possible to reduce this PVE and achieve better image quality by measuring and modeling the point spread function (PSF) and then accounting for it inside the reconstruction algorithm. In this work, we examined the response characteristics of the Metis^TM PET/CT system by acquiring ²²Na point source at different locations in the field of view (FOV) of the scanner and reconstructing with small pixel size for images to obtain their radial, tangential, and axial full-width half maximum (FWHM). An image-based model of the PSF model was then obtained by fitting asymmetric two-dimensional Gaussians on the ²²Na images. This PSF model determined by FWHM in three directions was integrated into a three-dimensional ordered subsets expectation maximization (3D-OSEM) algorithm based on a list-mode format to form a new PSF-OSEM algorithm. We used both algorithms to reconstruct point source, Derenzo phantom, and mouse PET images and performed qualitative and quantitative analyses. In the point source study, the PSF-OSEM algorithm reduced the FWHM of the point source PET image in three directions to about 0.67 mm, and in the phantom study, the PET image reconstructed by the PSF-OSEM algorithm had better visual effects. At the same time, the quantitative analysis results of the Derenzo phantom were better than the original 3D-OSEM algorithm. In the mouse experiment, the results of qualitative and quantitative analyses showed that the imaging quality of PSF-OSEM algorithm was better than that of 3D-OSEM algorithm. Our results show that adding the PSF model to the 3D-OSEM algorithm in the Metis^TM PET/CT system helps to improve the resolution of the image and satisfy the qualitative and quantitative analysis criteria.

1. Introduction

Positron emission tomography (PET) is a nuclear medical imaging technique and provides important information for disease diagnosis, therapeutic effect assessment, and new drug development [1]. Performing high-resolution PET imaging in small animals presents the challenge of achieving the same performance with target volumes that are orders of magnitude smaller than humans while maintaining adequate sensitivity and spatial resolution (1 mm or less) [2]. Despite significant advancements in the hardware and software of PET imaging systems, the spatial resolution of PET images is significantly lower than that of CT or MRI images [3]. There are many factors affecting the spatial resolution of PET, such as the size of the detector, the photon non-colinearity, the positron range and inter-crystal penetration [4,5,6]. These factors can lead to cross-contamination between voxels, referred to as the partial volume effect. In Ref. [7], improving the quality of the image by partial volume correction (PVC) was achieved. However, this method is not perfect.

The PET scanner consists of some rings of detectors arranged in a circular pattern. Due to this geometrical feature, distortions are introduced into the process of detection. If the photon enters the crystal from the center of the field of view (FOV), the line of response (LOR) is likely to be localized correctly. However, if the photon does not enter the crystal from the center of the FOV, there is a great possibility that the LOR will be positioned incorrectly. The further the photon is from the center of the FOV, the more likely it is that the LOR will be calculated incorrectly. This is because a photon travels through a crystal at an angle and then continues to travel through another crystal until it reaches its destination, and the result of this is an increase in distortion at points far from the center of the field of view [8]. A PSF modeling algorithm compensates for this distortion by incorporating measured PSFs at several million points in the FOV into its estimation step [9]. Recently, the information provided by the PSF and time-of-flight (TOF) has been expected to improve the spatial resolution and signal-to-noise ratio (SNR), respectively, of PET images [10,11]. In the past period of time, many researchers have become interested in the PSF model [12,13] and have received support and affirmation from other researchers and suppliers.

According to the cylindrical geometry of the PET detector arrangement, the full-width half maximum (FWHM) is reported along three orthogonal axes (radial, tangential, and axial). A particular point worth mentioning is that there is a spatial variation in the PSF across the FOV. The main reason for this is the depth-of-interaction (DOI) effect within the detector element [14]. There is a misclassification of photons traveling along lines that are more tangential to the FOV into a more central overlap line that extends throughout the FOV. This is because they enter the crystal surface of a particular detector element at a smaller angle, allowing them to penetrate neighboring crystals before they are detected. As a result, the PSF on the radial side of the FOV increases, while the PSF on the tangential side remains unchanged. All of these cases require that point sources be measured at different locations within the FOV to estimate the PSF [15].

Incorporating PET system response dates back to the introduction of filtered back-projection (FBP) algorithms [16,17]. A PSF-modeled system matrix contains more LORs, as each voxel corresponds to more measurement locations than a non-PSF-modeled system [18]. The PET technique is prone to partial volume effects (PVE), or image blurring, as well as signal spillover between adjacent functional regions [19]. Due to inconsistencies between the approximate (forward) projections of an image, an oversimplified projector converts the image signal into noise. In order to achieve high-quality reconstructions of images, it is necessary to develop an accurate model of the relationship between the image and projection space. This study aims to develop a PSF-OSEM reconstruction algorithm based on the combination of 3D-ordered subset expectation maximization (3D-OSEM) in list-mode format and PSF model to improve the spatial resolution and image quality of the associated PET system. To calculate the integral for each voxel, the continuous PSF model was discretized [20]. An effective measure of the spatial resolution of a scanner can be obtained by placing a point source in the scanner and acquiring scan data at varying locations along both the radial and axial directions [21,22]. As a part of this study, three sets of pre-experiments were designed using the 3D-OSEM algorithm in which a different number of iterations and subsets were selected to reconstruct ²²Na point sources. Then three different PSF models were obtained, and the accuracy of these three PSF models was compared to determine which PSF model was the most accurate for integration into the 3D-OSEM algorithm.

This paper presents the measurement and implementation of an approximate image model based on the response function of a small animal PET/CT scanner designed for the 3D-OSEM reconstruction algorithm in the Metis^TM PET/CT system developed by Shandong Madic Technology Co. Ltd. in China. We present measurements of the PSF obtained by acquiring and reconstructing ²²Na point sources. The response characteristics for this PET system were obtained by scanning ²²Na point sources at different locations in the scanner FOV and reconstructing them using the 3D-OSEM algorithm with five iterations and seven subsets. We analyzed the reconstructed images using the image analysis tool PMOD to calculate the radial, tangential, and axial FWHM of ²²Na point sources at different locations in the FOV, i.e., the parameters included in the PSF. By choosing the appropriate fitting function to fit the radial, tangential, and axial FWHM of all ²²Na point sources in scanned positions, we can determine the radial, tangential, and axial FWHM of all voxel points on the radial central slice in space. Taking into account the geometric characteristics of the PET system, we obtained the FWHM of all voxel points in scanner space in three directions by rotation. The value of the PSF of each voxel point was calculated by the PSF formula and stored in a three-dimensional matrix. It is this matrix that represents the PSF model, a three-dimensional asymmetric Gaussian distribution with three parameters. Lastly, the PSF model was incorporated into the 3D-OSEM algorithm based on the list-mode format to form the PSF-OSEM algorithm. This algorithm was evaluated on point source data, phantom data, and mouse images. The common 3D-OSEM algorithm (without PSF modeling) and the PSF-OSEM algorithm (with PSF modeling) were compared for visual inspection and quantitative analysis.

2. Materials and Methods

2.1. Materials

2.1.1. Point Source

In this paper, we used a point source of ²²Na with an activity of 1.22 MBq to measure the radial, tangential, and axial FWHM at different locations in the FOV. ²²Na is encapsulated in an acrylic cube with a side length of 1 cm and a diameter of 0.25 mm.

2.1.2. PET/CT System

This study was conducted using an integrated PET/CT system based on lutetium silicate scintillation crystal (LYSO), a device developed by Shandong Madic Technology Co. Ltd. in Shandong, China specifically for rodent imaging. An illustration of the appearance of the Metis^TM PET/CT system can be found in Figure 1a. This system has a cylindrical PET tomography scanner with 32 detector blocks arranged in four rings. The diameter of each ring is 129 mm, while the diameter of the effective ring is 81 mm. The axial range of the system is 122 mm. Each detector is equipped with 1152 LYSO crystals with a size of 0.943 × 0.943 × 10 mm³. The center-to-center distance between crystals is 1.028 mm. The crystal blocks are separated by a gap of 1.749 mm. A SiPM plate is attached to each crystal block. In Figure 1b, a simulated Metis^TM PET/CT device system is shown using GATE, a high-energy physics simulation toolkit based on Monte Carlo simulations.

2.2. Methods

2.2.1. PSF Measurement

The use of iterative reconstruction algorithms based on statistical models can maximize resolution recovery by accurately modeling the response of the system [23]. There are many LORs collected for each dataset, which poses formidable computational challenges for 3D iterative reconstruction. However, it is feasible to consider the response function of the detector on top of traditional iterative reconstruction methods such as OSEM. Due to multiple factors, including a large number of LORs and the difficulty of localizing the point source, it is unrealistic to conduct an ideal experiment using a collimated point source that emits photons along a single LOR to determine the degree of ambiguity in each LOR. To determine the system response matrix, PSF can be generated analytically, experimentally, and using Monte Carlo methods [24]. This study used experimental measurement methods.

2.2.2. PSF Estimation

Image spatial resolution deterioration is primarily caused by PVE, which can be mitigated by modeling the data for distortion before or during reconstruction. The PSF model is constructed during the reconstruction of the image in this experiment. It is always important to accurately model the PSF of a scanner during reconstruction to maximize its potential performance [23]. The voxels of the central slice are the only ones that require detailed modeling in tomography that have shift-invariant axial symmetry [2].

In this experiment, the radial, tangential, and axial FWHMs of the radial center slice were obtained by curve fitting the information of ²²Na point sources at different locations. By rotating the entire radial center slice, we can obtain the FWHM in three directions for all voxel points in the entire space.

In the PSF model, the three directions of FWHM in FOV correspond to the three spread parameters (radial spread parameter, tangential spread parameter, and axial spread parameter). In the field of optical imaging, if the PSF function is determined, it can be convolved with the acquired image. However, in practice, it is almost impossible to obtain the PSF accurately, and it can only be estimated through some experiments. In [20], the PSF function is defined as a Gaussian function, so the formula for obtaining the PSF value of this experiment is as follows:

$$PSF\left(x,y,z\right)=\frac{1}{{\left(2\pi \right)}^{\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}{\sigma }_r{\sigma }_t{\sigma }_a}{e}^{-\frac{x^2}{2{\sigma }_r^2}}{e}^{-\frac{y^2}{2{\sigma }_t^2}}{e}^{-\frac{z^2}{2{\sigma }_a^2}}$$

(1)

The coordinates (x, y, z) indicate the distance from each point to the center of the FOV in radial, tangential, and axial directions.

For the 3D PET images of ²²Na point sources collected at each location in space, we selected three 2D orthogonal planes (radial and tangential orthogonal, radial and axial orthogonal, and tangential and axial orthogonal) with maximum intensity voxels from each 3D image. Generally, the strongest voxel point is located at the center of the point source. We calculated the three spread parameter values for each ²²Na point source location by PMOD after confirming the three 2D orthogonal planes of the maximum intensity voxel. Each spread parameter was calculated twice. Finally, each spread parameter value was calculated as the average of the corresponding two values. Following this, spread parameters are estimated for each voxel in the space using an appropriate function. We selected the following fitting function based on the characteristics of the parameters:

$$f\left(x,z\right)=A_0{x}^2+A_1{z}^2+A_{2}xz+A_{3}x+A_{4}z+A_5$$

(2)

where x indicates the radial distance from the center of the FOV and z indicates the axial distance from the center of the FOV.

2.3. Experiments

2.3.1. Pre-Experiment

In order to measure the PSF value of each point in the space, we used the experimental measurement method. First, raw point source data were corrected for detector efficiency and dead time, and then the point source of ²²Na was placed at a fixed distance from the center of the FOV, and 12 axial positions were acquired (the point source to the center of the cross-axis FOV along the axial horizontal line at distances of −55 mm, −45 mm, −35 mm, −25 mm, −15 mm, −5 mm, 5 mm, 15 mm, 25 mm, 35 mm, 45 mm, 55 mm). In order to perform cross-axis FOV radial sampling along the horizontal line, tangential distances between the point source and the fixed distance are −30 mm, −25 mm, −20 mm, −15 mm, −10 mm, −5 mm, 0 mm, 5 mm, 10 mm, 15 mm, 20 mm, 25 mm, and 30 mm, in that order, after the axial position is fixed. For each ²²Na point source in the FOV, Figure 2 presents a model plot of the acquisition locations. The total number of points is 156.

2.3.2. Curve Fitting

By using the image processing and analysis software PMOD, we calculated the FWHM of the reconstructed data for each position in radial, tangential, and axial directions. According to the information on these three spread parameters obtained from the pre-experiment, we fit these data according to Equation (2). Figure 3 shows the FWHM fitting results of ²²Na point sources obtained at 156 locations. The FWHM and the fitting function can be matched well in all three directions at all positions except for a few special points.

By fitting the function, the values of the coefficients A₀, A₁, A₂, A₃, A₄, and A₅ can be obtained. Therefore, the three diffusion parameters can be calculated for each voxel point in the FOV of the scanner. Equation (1) can be used to calculate the PSF value for each voxel point in the whole space by adding the values of the parameters in the three directions. The PSF values were then stored in a 3D matrix (N × N × M), where N is the number of pixels in the radial or tangential direction of the reconstructed image, and M is the number of pixels in the axial direction of the reconstructed image.

2.4. Algorithm Optimization

The most widely used iterative algorithms in PET are the ML-EM (maximum-likelihood expectation maximization) algorithm and its accelerated version OSEM [25]. Statistical image reconstruction methods, such as OSEM, take into account the stochastic nature of the imaging process. These iterative algorithms have been shown to offer a better trade-off between noise and resolution in comparison to FBP [26,27]. The OSEM algorithm groups the projection data into an ordered sequence of subsets, and the standard EM algorithm is applied with a single pass through all subsets. This technique increases the convergence speed by a factor roughly equal to the number of subsets employed [26]. List mode is an efficient format for processing sparse datasets, such as dynamic or low-count studies. The list-mode 3D-OSEM algorithm formula can be described as follows [28].

$${\lambda }_j^{m,l}=\frac{{\lambda }_j^{m.l-1}}{N_j}{\displaystyle \sum }_{kϵS_l}{p}_{i_{k}j}\frac{A}{s_{i_k}+r_{i_k}+{{\displaystyle \sum }}_{b=1}^J{p}_{i_{k}b}{\lambda }_b^{m,l-1}}$$

(3)

where ${\lambda }_j^{m,l}$ indicates the 3D reconstructed image after $m$ iterations and $l$ subsets, and the subscript index of the voxel is $j=1,2\dots \dots J$; $s_l$ indicates that the event is divided into $s_l$ subsets, and the subscript $l$ is the number of subsets. For list mode, the subsets are formed according to the arrival time of the event [28]. $p_{ij}$ is the probability that the $j$th voxel produces a pair of gamma rays on the $i$th LOR; $s_{i_k}$ and $r_{i_k}$ are the scattering coefficient and the random coincidence coefficient of the LOR, where the $i_k$th event is located, respectively, and $A$ is the correction factor.

An accurate spatially variant resolution model for PET has been shown to reduce quantitative errors [29,30] and improve resolution by deconvolving the blurring [28]. This algorithm is modified by integrating the PSF model obtained from the PET system into the 3D-OSEM algorithm. Equation (3) is modified to Equation (4),

$${\lambda }_j^{m,l}=\frac{{\lambda }_j^{m.l-1}}{N_j}{\displaystyle \sum }_{kϵS_l}{p}_{i_{k}j}\frac{Ap{s}_{i_k}}{s_{i_k}+r_{i_k}+{{\displaystyle \sum }}_{b=1}^J{p}_{i_{k}b}\left[(·)\ast PSF\right]}$$

(4)

where $p{s}_{i_k}$ is the PSF value of voxel $b$ and * indicates the convolution. Each voxel contains contributions from collinear voxels. The size of the convolution kernel can be set according to the requirements. The convolution is calculated as shown below.

$${\left[(·)\ast PSF\right]}_{i,j,k}={{\displaystyle \sum }}_{r=-R}^R{{\displaystyle \sum }}_{t=-T}^T{{\displaystyle \sum }}_{a=-A}^A{\lambda }_{i+r}{\lambda }_{j+t}{\lambda }_{k+a}PSF\left[i+r,j+t,k+a\right]$$

(5)

We validated the proposed algorithm PSF-OSEM using point source, phantom data, and images of small animals. Additionally, a comparison of spatial resolution was conducted.

3. Validation

3.1. Point Source

Initially in the experiments, the point source of ²²Na with an activity of 1.22 MBq was acquired and reconstructed (3D-OSEM algorithm) at different radial and axial locations in the FOV in order to test the spatial resolution recovery introduced by PSF compensation. After integrating the PSF model into 3D-OSEM to form the PSF-OSEM algorithm, the same collected ²²Na point source data were reconstructed by the 3D-OSEM algorithm and the PSF-OSEM algorithm. In the experiment, we designed three different precision PSF models. The three models with different accuracy are combined with 3D-OSEM, respectively, to compare the quality of reconstructed images. According to the comparison results, the most accurate model was selected to be combined with 3D-OSEM and used in other experiments.

The number of iterations was based on a previous observation that the inclusion of resolution modeling with 3D-OSEM-PSF decreases the convergence rate of the algorithm [31]. Therefore, in order to compare the spatial resolution of these two algorithms, five iterations and seven subsets were set when reconstructing with the 3D-OSEM algorithm, and ten iterations and seven subsets were set when reconstructing with the PSF-OSEM algorithm. In addition, the size of the PSF kernel also has a significant impact on the reconstruction results, so we also compared the 3D-OSEM algorithm with the PSF-OSEM algorithm (with kernel sizes of three and five), and all three groups of algorithms were performed for five iterations and seven subsets.

3.2. Derenzo Phantom

Image quality was evaluated with the Derenzo phantom [32]. Figure 4 [32] illustrates six groups of aperture (0.5, 0.6, 0.7, 0.8, 0.9, 1.0 mm) of the Derenzo phantom. The center-to-center distance between adjacent rods was twice the diameter of the rod in the same group. PET/CT with fluorine-18 (¹⁸F)-fluorodeoxyglucose (FDG) is useful for the detection and staging of various malignant tumors, monitoring of their response to therapy, and prognostic stratification [33,34,35,36,37,38]. All the rods of the Derenzo phantom were filled with ¹⁸F-FDG mixed with saline. The injected dose was 3.71 MBq. Then the phantom was placed in the fourth ring of the scanner (as near as possible to the end of the PET axial FOV) and scanned for 20 min.

3.3. Small Animal Studies

Animal studies were approved by the Laboratory Animal Ethics Committee of Xuzhou Medical University (Process number for animal experiments: 201706w010). The solution of ¹⁸F-FDG with an activity of 18.95 MBq was injected into a healthy mouse through the tail vein. The weight of the mouse was 43.78 g, its length was 11 cm, and the mouse was placed in the center of the axial FOV. During the experiment, the mouse was anesthetized with a mixture of oxygen and isoflurane. During the scanning, the laser was applied to the tail of the mouse, and multi-bed scanning was conducted for 30 min.

3.4. Reconstruction Parameters

The collected point source and the Derenzo phantom data were reconstructed by the 3D-OSEM algorithm and PSF-OSEM algorithm with two different sets of parameters.

We selected different parameters for the reconstructions in order to provide more data to support the conclusion that PSF-OSEM can improve the quality of the images and that PSF modeling can slow down the iterative convergence process. In PSF-OSEM, several PSF kernel sizes are available (in this case, three and five kernels). The following are the parameters that are set in the reconstructions:

• Point sources. For 3D-OSEM algorithm reconstruction, the 3D image matrix was 257 × 257 × 389, subsets were 7, the iteration was 5, the reconstruction FOV diameter was 81 mm, the pixel size was 0.314 mm, and the Gaussian post-filter was FWHM = 2 mm. For PSF-OSEM algorithm reconstruction, the 3D image matrix was 257 × 257 × 389, the subsets were 7, the iteration was 10, the reconstruction FOV diameter was 81 mm, the pixel size was 0.314 mm, and the Gaussian post-filter was FWHM = 2 mm;
• Derenzo phantom. Reconstruction parameters are the same for both algorithms: the 3D image matrix was 257 × 257 × 389; the iterations were 2, 5, 7, 10, 12, 15, 20; the subsets were 5, the reconstructed FOV diameter was 81 mm; pixel size was 0.314 mm; and Gaussian post-filter was FWHM = 2 mm;
• Small animal studies. Reconstruction parameters were the same for both algorithms: the 3D image matrix was 257 × 257 × 389, the subsets were 5, the iterations were 7, the reconstructed FOV diameter was 81 mm, the pixel size was 0.314 mm, and the Gaussian post-filter was FWHM = 2 mm.

3.5. Data Analysis

The purpose of the data analysis is to evaluate whether the improved PSF-OSEM algorithm can improve the resolution of the image.

Point source. Two different reconstruction algorithms were compared with respect to the one-dimensional distribution of FWHM along the radial, tangential, and axial directions. Furthermore, the influence of the kernel size in the PSF-OSEM algorithm on the FWHM was also investigated.

Derenzo phantom. We evaluated image quality by visual assessment first and then quantitatively analyzed the images based on three parameters: contrast-to-noise ratio (CNR), contrast, and contrast recovery (CR_hot). We regarded the smallest hot rod that clearly distinguished the PET image of the Derenzo phantom as a criterion for visual inspection. We outlined a circular region of interest (ROI) with a diameter of 0.8 mm on a 0.8 mm hot bar and then a background region with a diameter of 1.6 mm. In each case, the 133rd to 138th layers were selected. The three evaluation metrics are formulated as follows.

• CNR:

  $$CNR=\frac{S_{mean}-B{D}_{mean}}{S{D}_{background}}$$

  (6)
• Contrast:

  $$Contrast=\frac{S_{mean}}{B{D}_{mean}}$$

  (7)
• CR_hot,:

  $$C{R}_{hot}=\frac{\raisebox{1ex}{$S_{mean}$}\!\left/ \!\raisebox{-1ex}{$B{D}_{mean}$}\right.-1}{R-1}$$

  (8)

  where $S_{mean}$ is the mean signal intensity of the ROI in the Derenzo phantom, $B{D}_{mean}$ is the mean of the intensity of the background region, and $S{D}_{background}$ is the standard deviation of the background. $R$ is the real signal-to-background ratio.

Small animal studies. Qualitative and quantitative analyses were performed on PET images of the mouse reconstructed by the two algorithms with seven iterations and five subsets.

4. Results

4.1. ²²Na Point Source

Figure 5 illustrates the one-dimensional representation of the FWHM along the radial, tangential, and axial axes of the ²²Na point source images reconstructed by the two algorithms at a distance of 5 mm and 15 mm axially from the center of the FOV, respectively. In three different directions, the FWHM of the PSF-OSEM algorithm (the PSF convolution kernel here is three) is significantly lower than the FWHM of the original 3D-OSEM algorithm. As the radial distance increases, the radial FWHM of the two algorithms also increases. However, tangential and axial FWHM is not significantly related to radial distance. In Figure 6, FWHM is presented in three dimensions after reconstruction of the same point source data at a distance of 35 mm axially from the center of the FOV using the 3D-OSEM algorithm and the PSF-OSEM algorithm with two different PSF kernel sizes (three and five). The reconstruction parameters of the three reconstruction algorithms are all five iterations and seven subsets. It should be noted that the PSF-OSEM algorithm has a more uniform spatial resolution over the entire FOV, and the size of the PSF kernel was found to have a significant impact on the resolution. Additionally, it was found that the PSF-OSEM algorithm performs significantly better with a PSF kernel of five than with a PSF kernel of three. Therefore, it appears that the PSF-OSEM algorithm is useful in improving the spatial resolution of the PET/CT system, and when the kernel is five, the spatial resolution is better.

4.2. Derenzo Phantom

Compared to the PET images reconstructed by the 3D-OSEM algorithm, the PSF-OSEM algorithm improves the visualization of the PET images of the Derenzo phantom at the hot rod of 0.6 mm. Figure 7 shows PET images of the Derenzo phantom based on two different algorithms. For both algorithms, the reconstruction parameters are 5 subsets and 5 iterations, 7 iterations, 10 iterations, 12 iterations, 15 iterations, and 20 iterations. The same number of iterations and subsets are presented with different visual effects by different reconstruction algorithms. In Figure 7, it can be seen that the PET images reconstructed with the PSF-OSEM algorithm are clearer than those reconstructed with the 3D-OSEM algorithm. The hot rod of 0.6 mm in the image of the Derenzo phantom reconstructed with the PSF-OSEM algorithm can be clearly identified, while the 3D-OSEM algorithm does not distinguish the hot rod of 0.6 mm more easily than the PSF-OSEM algorithm. However, at the right edge of the Derenzo phantom, the noise of PET images reconstructed by the PSF-OSEM algorithm is higher than that of the 3D-OSEM algorithm. Visually, the 3D-OSEM algorithm reconstructs better image quality from 7 iterations and 5 subsets to 10 iterations and 5 subsets. The PSF-OSEM algorithm reconstructs higher-quality images from 10 iterations and 5 subsets to 15 iterations and 5 subsets. When the number of iterations increases, the images become blurred. Due to the different numbers of iterations, the metrics used to measure their image quality can be affected. According to Figure 8, there is a significant difference between the values of CNR, contrast, and CR_hot under different reconstruction algorithms. The PET images of the Derenzo phantom reconstructed by the PSF-OSEM algorithm have a higher CNR, contrast, and CR_hot than those reconstructed by the 3D-OSEM algorithm in general. The CNR and Contrast of the 3D-OSEM algorithm are only higher than those of the PSF-OSEM algorithm when the number of iterations is small.

4.3. Small Animal PET Images Analysis

Figure 9 shows the PET images of a healthy mouse reconstructed with the 3D-OSEM algorithm and PSF-OSEM algorithm, respectively. Figure 9b shows that the mouse’s brain, heart, bladder, and other organizational structures have higher contrast than those in Figure 9a, i.e., the structures and their “functional” morphology are more pronounced than those shown in Figure 9a. According to the quantitative analysis of the Derenzo phantom, when the number of iterations is greater than 7, its three analysis indicators show that the reconstructed image quality of the PSF-OSEM algorithm is better than that of the 3D-OSEM algorithm. We drew the ROI at the heart and bladder and analyzed them with the same three indicators as the Derenzo phantom.

Table 1. Three quantitative indexes of PET mouse images reconstructed using 3D-OSEM and PSF-OSEM algorithms.

	3D-OSEM			PSF-OSEM
	CNR	Contrast	CRhot	CNR	Contrast	CRhot
Heart	3.97	7.406	1.784	4.045	7.642	1.818
Bladder	3.878	11.582	3.262	4.196	13.105	3.416

shows that the three indicators of the PET image reconstructed by the PSF-OSEM algorithm are greater than those of the 3D-OSEM algorithm.

5. Discussion

PET/CT hybrid systems are currently a major trend, and the combination of PET and CT allows for more comprehensive image information. Several improvements have been made to the hardware of PET systems in order to achieve high-resolution images. Despite this, the imaging effect of PET is still not ideal due to the influence of a variety of physical factors. By modeling an accurate PSF model, we can effectively address the problem of image quality degradation caused by physical factors. Based on theoretical considerations, PSF can be effective in improving spatial resolution and determining the position of the response lines in space. According to the literature, there are three main methods to obtain PSF models: analytical modeling [23,39,40], simulation [41], and experimental measurement [30,42,43]. In this paper, experimental measurements were made to obtain the estimation of PSF. For the purpose of this experiment, an imaging approach, i.e., a PSF model coupled with a 3D-OSEM algorithm, was utilized to aid in the reconstruction of the images. PSF estimation is dependent on the reconstruction algorithm used when considering PSF at the image level. The advantage of this approach is that it is not only simple, but it also allows for improved accuracy and a better match between pixel size and PSF size. Moreover, it is more appropriate to estimate PSF based on image levels when dealing with list-mode data formats [20].

The purpose of this study is to propose a PSF-OSEM algorithm based on the 3D-OSEM algorithm in the Metis^TM PEC/CT system developed by Shandong Madic Technology Co., Ltd. in Shandong, China. PSF values for each voxel are based on three spread parameters: radial spread parameter, tangential spread parameter, and axial spread parameter. For this system, it is essential that the response function be measured carefully, a small point source should be selected, and relevant factors should be considered in order to obtain an accurate response. Among these factors is the size of the source, its type of isotope, the medium surrounding it (air, water, warm radioactive background), and the number and location of measurements so that all possible responses of the system can be considered with respect to the detector configuration [30,44]. Because the point source was located in the air, the convergence of 3D-OSEM was fast [31]. To measure the PSF values of this system, this experiment was performed using a ²²Na point source with an activity of 1.22 MBq, a sphere with a diameter of 0.25 mm, encapsulated in an acrylic cube with a side length of 1 cm. In many ways, the characteristics of this type of point source are similar to those of the ¹⁸F point source in water. One hundred and fifty-six point source measurements were sampled uniformly in the FOV of the scanner, which was a very important component of the study. We selected five iterations and seven subsets for image reconstruction, and a Gaussian filter with an FWHM of 2 mm was used as a post-smoothing filter. As a result of each reconstructed image, three orthogonal two-dimensional planes (radial and tangential, radial and axial, and tangential and axial) with maximum intensity voxels were extracted, and each plane was fitted with an appropriate function in order to obtain diffusion parameters for all points in the space. We also designed different PSF models with varying numbers of iterations and subsets during the pre-experimental stage, but there was no significant difference in the improvement of PET image quality. This pre-experiment is not reflected in this paper.

It was validated on point source data, phantom models, and mouse images to verify whether the new algorithm could function correctly and improve the image quality of PET by incorporating PSF information into the anterior-posterior projection of the original 3D-OSEM algorithm. By using the same acquisition dataset, the imaging quality of the two algorithms before and after modification was compared. Figure 5 shows that at the position near the FOV center of the scanner, the radial FWHM using the 3D-OSEM algorithm is between 0.9 mm and 1.0 mm, while the radial FWHM using the PSF-OSEM algorithm is between 0.6 mm and 0.8 mm. Outside the radial distance of the FOV center of 20 mm, the radial FWHM using the 3D-OSEM algorithm is between 1.5 and 3.1 mm, while the radial FWHM using the PSF-OSEM algorithm is between 0.9 and 2.8 mm. Additionally, the radial FWHM of the PSF-OSEM algorithm is higher than 1.6 mm only at 30 mm. Regardless of how the radial distance changes, the tangential and axial FWHM changes little. Tangential and axial FWHM using the PSF-OSEM algorithm is much lower than the 3D-OSEM algorithm. The tangential and axial FWHM using the 3D-OSEM algorithm is about 0.9 mm and 1.0 mm, while the tangential and axial FWHM using the PSF-OSEM algorithm are about 0.6 mm and 0.7 mm. The FWHM in three directions using the PSF-OSEM algorithm is smaller than that using the 3D-OSEM algorithm, and the size of the PSF kernel also affects the FWHM value in three directions. Figure 6 shows that when the PSF kernel = 5, the FWHM in three directions is the minimum, and the change in FWHM value is relatively stable. The PSF-OSEM algorithm (with a kernel of 5) produces a more uniform radial FWHM in the FOV compared to the 3D-OSEM algorithm. The FWHM of point sources reconstructed using the PSF-OSEM algorithm is significantly lower than the 3D-OSEM algorithm in three directions, and it is possible to overcorrect the PSF resolution when evaluating air point sources due to unacceptable noise characteristics, which makes it impractical for clinical application. A series of experiments were performed with the Derenzo phantom to validate the conclusion that adding PSF to the 3D-OSEM algorithm improves spatial resolution. As shown in Figure 7, in the qualitative analysis of Derenzo phantom, the PET image using the PSF-OSEM algorithm is clearer than the PET image using the 3D-OSEM algorithm. In the PET images of the Derenzo phantom reconstructed by the two algorithms, hot rods above 0.6 mm can be distinguished, but it is easier for the PSF-OSEM algorithm to distinguish the smallest visible hot rod (0.6 mm). As the number of iterations increases, the PET image reconstructed by the 3D-OSEM algorithm gradually blurs, while the PSF-OSEM gradually becomes clear, which is related to the PSF model slowing down the iterative convergence speed. It should be noted that the increased number of iterations implies that there is more noise. From the PET images of the Derenzo phantom, it can be seen that the parameters of a good PET image reconstructed by the 3D-OSEM algorithm are from 7 iterations and 5 subsets to 10 iterations and 5 subsets, while the parameters of a good PET image reconstructed by PSF-OSEM algorithm are from 10 iterations and 5 subsets to 15 iterations and 5 subsets. Next, quantitative analyses were made on the PET images of the Derenzo phantom reconstructed by two algorithms: CNR, contrast, and CR_hot. It is obvious from Figure 8 that the three indicators of PSF-OSEM are higher than those of 3D-OSEM. The CNR and contrast of the PSF-OSEM algorithm are only lower than that of the 3D-OSEM algorithm when the number of iterations is small. In the mouse experiment, PET imaging using the PSF-OSEM algorithm had higher contrast than that using the 3D-OSEM algorithm, and the results of the quantitative analysis show that the three indicators of the PET image reconstructed by the PSF-OSEM algorithm are higher. In the ROI of the heart, although the three index parameters of the reconstructed image using the PSF-OSEM algorithm are higher, the numerical difference is small, which is where we need to improve further in the later stage.

The PSF-OSEM algorithm improves image quality for PET images. However, the selection of reconstruction parameters also determines the quality of PET images. The next step is to select the optimal imaging parameters for the PSF-OSEM algorithm. At the same time, because of the PSF model, the time of image reconstruction becomes longer. Reducing reconstruction time is also the focus of the next study.

6. Conclusions

By measuring the ²²Na point source using the 3D-OSEM algorithm reconstruction and image analysis tool PMOD, the FWHM in three directions of 156 positions in the scanner FOV is obtained, and then the PSF model of the entire Metis^TM PET/CT system is established by fitting with appropriate functions. The modified PSF model is integrated into the 3D-OSEM algorithm to form the PSF-OSEM algorithm. Through experimental verification, the modified PSF-OSEM algorithm can improve the spatial resolution of the system and meet the qualitative and quantitative analysis indicators. The reconstruction time of this algorithm is longer than the original 3D-OSEM algorithm, but considering the improvement of image quality, the sacrifice of time is worth it.