A Time-of-Flight and GPU-Based Method for PET Scattering Correction and Image Reconstruction

Abstract

This invention discloses a PET scattering correction and image reconstruction method based on time-of-flight (TOF) and GPU, belonging to the field of medical imaging technology. This invention uses a TOF histogram image as the initial radiation source distribution for scattering estimation for scattering correction. The computational overhead of the TOF histogram image is extremely low, and the initial scattering estimation can be executed in parallel with the sensitivity image calculation process, thereby effectively improving the computational efficiency of the overall reconstruction process. This invention has strong flexibility; it can be modularly integrated into EM-type reconstruction algorithms without requiring significant modifications to existing reconstruction frameworks, facilitating implementation. An atomically driven multi-GPU parallel MC computing framework solves the load balancing problem in multi-GPU computing, fully leveraging the performance advantages of multi-GPU computing resources. While meeting clinical needs in terms of computational efficiency, it can accurately estimate the scattering distribution, including multiple scattering coincidence events. The method of this invention is simple and easy to understand. This invention is applied to PET scattering correction and image reconstruction.

Description

Technical Field

[0001] This invention relates to medical imaging technology, specifically to a PET scattering correction and image reconstruction method based on time of flight and GPU. Background Technology

[0002] Positron emission tomography (PET) is an important functional imaging technique widely used in oncology, neuroscience, and cardiovascular disease. During PET imaging, β... + The positrons produced by the decay of radioactive nuclides annihilate electrons in human tissue, emitting a pair of gamma photons with an energy of 511 keV and an approximately 180° directional orientation. PET detection systems generate coincidence events by detecting the coincidence of these photon pairs. Using appropriate image reconstruction algorithms, beta photons can be accurately and quantitatively reconstructed from a large number of coincidence events. + The spatial distribution of radionuclides within the body provides quantitative information on the distribution and metabolism of radiopharmaceuticals for clinical diagnosis and scientific research. However, some gamma photons undergo Compton scattering in human tissue or detectors before being detected, altering their energy and propagation direction, resulting in scatter coincidence events. Scatter coincidence events significantly reduce image contrast, increase background noise, and affect quantitative accuracy; therefore, scatter correction is an indispensable step in PET image reconstruction.

[0003] Currently, the most commonly used scattering correction method in mainstream PET equipment integrates image domain scattering estimation within the Expectation Maximization (EM) algorithm framework. Through joint iterative optimization, scattering artifacts in the reconstructed image are eliminated. Image domain scattering estimation methods include Monte Carlo (MC) simulation and Single-scatter simulation (SSS) methods. Among them, the MC method is considered the "gold standard" for scattering estimation due to its accurate physical modeling; however, its computational cost is enormous, and its computational efficiency is insufficient to meet clinical needs, limiting its application in practical systems. EM algorithms are the "gold standard" for PET image reconstruction, mainly divided into two categories: Maximum Likelihood Expectation Maximization (MLEM) and Ordered Subset Expectation Maximization (OSEM). OSEM is an optimized version that accelerates MLEM image reconstruction by partitioning matching events into subsets. Figure 1 shows a typical flow of an EM-type reconstruction algorithm that integrates scattering event estimation and correction functions. In an image reconstruction process involving N iterations, sensitivity image calculation, the initial EM reconstruction without scattering correction, and the nth (n=2,3,…,N, where N is the iteration number) EM reconstruction with scattering correction are executed sequentially, with data dependencies between each step. The activity distribution image without scattering correction obtained from the initial EM reconstruction is used as the initial radiation source distribution input to the scattering estimation program to roughly estimate the distribution of scattering events and calculate the corresponding scattering factor. This scattering factor is then input into the second EM reconstruction to achieve preliminary scattering correction. Subsequently, the activity distribution image with scattering correction output from each EM reconstruction is used as a new input to the scattering estimation program to finely update the scattering factor and perform the next EM reconstruction with scattering correction. This process iterates until a preset termination condition is met, and the activity distribution image output from the last iteration is the reconstruction result. With the development of long-axis field-of-view PET (LAFOV-PET) systems, the number of coincidence events acquired in a single scan has increased significantly, placing higher demands on the computational efficiency of EM-type image reconstruction algorithms. Simultaneously, the proportion of coincidence events from multiple scattering events has increased significantly, making it difficult for the SSS method based on the single scattering assumption to achieve accurate scattering estimation. Therefore, it is necessary to re-examine the overall process of EM-type image reconstruction algorithms to balance computational efficiency and scattering estimation accuracy.

[0004] In recent years, Time-of-Flight (TOF) detection technology has become a standard feature of commercial PET systems, with some advanced PET systems achieving TOF resolutions of around 200 ps. Existing reconstruction techniques incorporate TOF information into weighted or inverse projection or system matrix modeling during EM reconstruction, significantly improving the signal-to-noise ratio and quantitative accuracy of the reconstructed images. However, the potential of TOF information in constructing the initial radioactive source distribution for scatter estimation has not yet been fully explored. Current methods typically use the uncorrected activity distribution image obtained from the initial EM reconstruction as the initial radioactive source distribution for scatter estimation. This not only requires an additional EM reconstruction but also prevents the initial scatter estimation process from being parallelized with other steps in the reconstruction workflow, resulting in low overall computational efficiency. Summary of the Invention

[0005] To address the problems of existing technologies, this invention proposes a PET scattering correction and image reconstruction method based on time-of-flight and GPU. Compared to the traditional method of using the first EM reconstruction image without scattering correction as the initial radiation source distribution, the computational overhead of the TOF histogram image is significantly reduced, and the first scattering estimation can be performed in parallel with the sensitivity image calculation process, thereby effectively improving the computational efficiency of the overall reconstruction process. Furthermore, to achieve high-precision scattering estimation including multiple scattering coincidence events, this invention employs the MC method for scattering modeling. To further improve the computational efficiency of the MC method and meet clinical needs, this invention proposes an atom-driven multi-GPU (Graphic Processing Unit) parallel MC computing framework to solve the load balancing problem in multi-GPU computing and fully leverage the performance advantages of multi-GPU computing resources.

[0006] The PET scattering correction and image reconstruction method based on time of flight and GPU of the present invention includes the following steps:

[0007] 1) The PET device collects coincidence events, calculates TOF information based on the coincidence events, and calculates the most likely annihilation location based on the geometric information of the PET and the TOF information; calculates the weight of the coincidence event based on the normalized information of the detection efficiency and the attenuation image; and calculates the TOF histogram image based on the most likely annihilation location, coincidence event type and weight.

[0008] 2) The TOF histogram image is used as the initial radiation source distribution for scattering estimation; based on the attenuation image and the initial radiation source distribution, a preliminary scattering estimation of the scattering factor is performed to obtain the preliminary estimated scattering factor that matches the event; in the process of preliminary scattering estimation, atomic-driven multi-GPU parallel computing is used.

[0009] 3) Based on the geometric information of PET, the normalized information of detection efficiency, and the attenuation image, the sensitivity image is calculated to obtain the sensitivity image; Steps 1) and 2) are calculated sequentially, and Steps 1) and 2) are calculated in parallel and synchronously with Step 3).

[0010] 4) The EM algorithm is used for image reconstruction. The image reconstruction process includes all physical corrections. Among them, the scattering correction in the physical correction first calculates the activity distribution image with the initial scattering correction based on the initial estimated scattering factor, coincidence event, and sensitivity image. Then, the scattering factor is finely estimated based on the activity distribution image to obtain the finely estimated scattering factor. The finely estimated scattering factor is introduced into the subsequent iteration to update the activity distribution image, realizing the coordinated iterative update of the scattering factor and the activity distribution image. The reconstructed image is obtained after the iteration converges.

[0011] In step 1), the geometric information of the PET sensor represents the spatial distribution of the PET detector crystal. The attenuation image is calculated from the CT image, providing information on the attenuation of gamma photons by the human anatomical structure. The TOF histogram calculation includes the following steps:

[0012] a) Record matching events:

[0013] Each coincidence event includes: the crystal position of the fast photon. Fast Photon Timestamp Crystal position of slow photons Slow photon timestamp and matching event type The event types include instantaneous events and delayed events. For instantaneous events... =0, for delayed events =1;

[0014] b) Calculate TOF information:

[0015] Read the timestamps of two single-photon events from the matching events, and calculate the Time-of-Flight (TOF) information based on the timestamps of the two single-photon events. for:

[0016]

[0017] c) Calculate the most likely annihilation location:

[0018] Based on the TOF information and the crystal positions of the two single-photon events (i.e., the geometric information of the PET), the most likely annihilation position is calculated as follows:

[0019]

[0020] Where c represents the speed of light. The most likely annihilation location;

[0021] d) Calculate the weights:

[0022] Based on the normalized information of PET detection efficiency, the normalized correction factor corresponding to the matching event is calculated as follows: Perform orthographic projection calculations on the attenuation image to determine the attenuation correction factor corresponding to the event. Then calculate the weight corresponding to the matching event. for:

[0023]

[0024] e) Calculate voxel values:

[0025] Given a TOF histogram image X, where all voxels are initially set to 0, update the voxel values ​​of the TOF histogram image X using all matching events; if the most probable annihilation location... The first point that falls on the TOF histogram image X In individual elements, then according to the event... and weight Update # Individual physiology value :

[0026]

[0027] The calculation continues until all coincident events are calculated, resulting in a TOF histogram image. This calculation method can also be efficiently implemented in parallel on a GPU, where j=1,..,J, and J is the number of voxels in the TOF histogram image X.

[0028] In step 2), the TOF histogram image is used as the initial radiation source distribution for scattering estimation. Compared with the traditional method of obtaining the initial radiation source distribution through EM reconstruction, this significantly reduces the computational cost. The MC method or SSS method is used to perform a preliminary scattering estimation of the scattering factors. Given the poor accuracy of the SSS method, this invention employs an atom-driven multi-GPU parallel MC computing framework to perform scattering estimation. While maintaining the accuracy of the MC method's physical modeling, it fully utilizes the powerful parallel computing capabilities of multiple GPUs to achieve good load balancing, enabling high-precision scattering distribution estimation to be completed within a clinically acceptable timeframe.

[0029] In multi-GPU parallel computing, if the total number of radioactive atoms allocated to each GPU and their spatial distribution are inconsistent, it can easily lead to unbalanced computational load, thereby reducing overall computational efficiency and wasting computational resources. To address this, this invention proposes an atom-driven multi-GPU parallel MC computing framework. This framework effectively solves the load imbalance problem by allocating tasks at the atomic level with fine granularity, fully leveraging the parallel advantages of multi-GPU computing resources to achieve efficient and accurate scattering estimation.

[0030] The atomic-driven multi-GPU parallel MC computing framework includes the following steps:

[0031] a) Convert the TOF histogram image X into a radioactive atom distribution image Y. The conversion relationship is as follows:

[0032]

[0033] in, and The first two images represent the TOF histogram image X and the radioactive atom distribution image Y, respectively. Individual phenotypic value The total number of radioactive atoms to be simulated;

[0034] b) Let the PET acquisition time simulated by the MC method be... The total number of radioactive atomic decays simulated by the MC method is M:

[0035]

[0036] in, The half-life of a radioactive nuclide; and The value needs to be set by the user; generally, for scattering estimation, in order to make the MC method have a sufficiently low statistical variance, the number of radioactive atomic decays M calculated according to the above two formulas should be within 10. 8 Order of magnitude;

[0037] c) Employing an atomic-level load balancing strategy, by traversing all voxels of the radioactive atom distribution image Y, the set of radioactive atoms allocated to each GPU is computed in parallel across multiple GPUs with a total of K GPUs. The allocation relationship is as follows:

[0038]

[0039] in, This represents the set of radioactive atoms allocated to each GPU. In the The value at the individual element, This indicates a rounding down operation; the distribution is carried out according to the above formula to ensure that the radioactive atoms obtained by each GPU are consistent in terms of the total number and spatial distribution, so that the computing load undertaken by each GPU is basically the same, thereby achieving load balancing in the multi-GPU computing process and improving the overall efficiency of multi-GPU parallel computing.

[0040] d) Each GPU performs a preliminary scattering estimate based on the attenuation image according to the assigned set of radioactive atoms, and obtains a preliminary estimated scattering factor that matches the event.

[0041] In step 3), the sensitivity image is calculated using the standard procedure of the EM class algorithm.

[0042] In step 4), each iteration yields a scattering-corrected activity distribution image. This image is then used to perform a detailed scattering estimation of the scattering factor, resulting in a finer scattering factor. This is used for the next iteration. During the iterations, the scattering factor and the activity distribution image become increasingly accurate. The activity distribution image obtained after the iterations converge is the reconstructed image. The detailed scattering estimation employs a multi-GPU parallel McAfee (MC) computing framework. After N iterations, the iterations converge. The number of iterations N depends on the specific task and is typically N < 10.

[0043] This invention breaks the strict data dependency between sensitivity image calculation, the initial EM reconstruction without scattering correction, and subsequent EM reconstructions with scattering correction. This allows the TOF histogram image calculation and preliminary scattering estimation steps to be executed in parallel with sensitivity image calculation, thereby effectively improving the overall computational efficiency of the algorithm. It can save the computation time of one EM reconstruction without scattering correction and one scattering estimation.

[0044] Advantages of this invention:

[0045] This invention presents a scattering correction and image reconstruction method that uses Time-of-Flight (TOF) histogram images as the initial radiation source distribution for scattering estimation. The computational cost of TOF histogram images is extremely low, and the preliminary scattering estimation can be performed in parallel with the sensitivity image calculation process, thereby effectively improving the overall computational efficiency of the reconstruction workflow. This invention is applicable not only to scattering estimation using the Monte Carlo (MC) method but also to scattering estimation using the Saxophone Segmentation (SSS) method, demonstrating strong flexibility. It can be modularly integrated into EM-type reconstruction algorithms without requiring significant modifications to existing reconstruction frameworks, facilitating implementation. An atom-driven multi-GPU parallel MC computation framework solves the load balancing problem in multi-GPU computation, fully leveraging the performance advantages of multi-GPU computing resources. While meeting clinical needs in terms of computational efficiency, it can estimate the scattering distribution, including multiple scattering coincidence events, with high accuracy. The method is simple and easy to understand. This invention is applied to PET scattering correction and image reconstruction. Attached Figure Description

[0046] Figure 1 A typical flowchart of an existing EM-type reconstruction algorithm that includes scattering correction;

[0047] Figure 2 This is a flowchart of an embodiment of the time-of-flight and GPU-based PET scattering correction and image reconstruction method of the present invention;

[0048] Figure 3 This is a schematic diagram illustrating the histogram image calculation principle of the time-of-flight and GPU-based PET scattering correction and image reconstruction method of the present invention;

[0049] Figure 4 The TOF histogram is obtained by the PET scattering correction and image reconstruction method based on time of flight and GPU of the present invention.

[0050] Figure 5 The flowchart shows the multi-GPU parallel MC computing framework of the PET scattering correction and image reconstruction method based on time of flight and GPU of the present invention.

[0051] Figure 6 Images reconstructed using conventional methods and the method of this invention. Detailed Implementation

[0052] The present invention will be further described below with reference to the accompanying drawings and specific embodiments.

[0053] This embodiment of the present invention provides a time-of-flight and GPU-based PET scattering correction and image reconstruction method, such as... Figure 2 As shown, it includes the following steps:

[0054] 1) PET equipment acquires coincidence events. In this embodiment, the NEMA Image Quality phantom (hereinafter referred to as NEMA IQ phantom) defined in the NEMA NU 2-2018 standard is used. PET imaging simulation is performed using GATE software, with the TOF resolution set to 250 ps during simulation. The GATE software outputs a list of coincidence events, totaling 1,079,001,233 coincidence events, including transient and delayed events. The geometric information of the PET represents the spatial distribution of the PET detector crystal. The attenuation image is calculated from the CT image, providing information on the attenuation of γ-photons by the human anatomical structure. The TOF histogram is calculated as follows: Figure 3 As shown:

[0055] a) Record matching events:

[0056] Each coincidence event includes: the crystal position of the fast photon. Fast Photon Timestamp Crystal position of slow photons Slow photon timestamp and matching event type This includes instantaneous events and delayed events. For instantaneous events... =0, for delayed events =1;

[0057] b) Calculate TOF information:

[0058] Read the timestamps of two single-photon events from the matching events, and calculate the Time-of-Flight (TOF) information based on the timestamps of the two single-photon events. for:

[0059]

[0060] c) Calculate the most likely annihilation location:

[0061] Based on the TOF information and the geometric information of the crystal positions of the two single-photon events (i.e., PET), the most likely annihilation location is calculated as follows:

[0062]

[0063] Where c represents the speed of light. The most likely annihilation location;

[0064] d) Calculate the weights:

[0065] Based on the normalized information of PET detection efficiency, the normalized correction factor corresponding to the matching event is calculated as follows: Perform orthographic projection calculations on the attenuation image to determine the attenuation correction factor corresponding to the event. Then calculate the weight corresponding to the matching event. for:

[0066]

[0067] e) Calculate voxel values:

[0068] Given a TOF histogram image X, where all voxels are initially set to 0, update the voxel values ​​of the TOF histogram image X using all matching events; if the most probable annihilation location... The first point that falls on the TOF histogram image X In individual elements, then according to the event... and weight Update # Individual physiology value :

[0069]

[0070] This continues until all matching events have been calculated, resulting in a TOF histogram image, such as... Figure 4 As shown, this calculation method can also be efficiently implemented in parallel on a GPU, j=1,..,J, where J is the number of voxels in the TOF histogram image X;

[0071] 2) The TOF histogram image is used as the initial radiation source distribution for scattering estimation; based on the attenuation image and the initial radiation source distribution, the MC method is used to perform a preliminary scattering estimation of the scattering factor, obtaining a preliminary estimated scattering factor that matches the coincidence event; in the process of preliminary scattering estimation, an atom-driven multi-GPU parallel MC computing framework is used, such as... Figure 5 As shown:

[0072] a) Given that the MC method uses the decay process of radioactive atoms as the basic simulation unit, the TOF histogram image X is first converted into a radioactive atom distribution image Y. The conversion relationship is as follows:

[0073]

[0074] in, and The first two images represent the TOF histogram image X and the radioactive atom distribution image Y, respectively. Individual phenotypic value The total number of radioactive atoms to be simulated;

[0075] b) Let the PET acquisition time simulated by the MC method be... The total number of radioactive atomic decays simulated by the MC method is M:

[0076]

[0077] in, The half-life of a radioactive nuclide; and The value needs to be set by the user; generally, for scattering estimation, in order to make the MC method have a sufficiently low statistical variance, the number of radioactive atomic decays M calculated according to the above two formulas should be within 10. 8 Order of magnitude;

[0078] c) An atomic-level load balancing strategy is adopted. By traversing all voxels of the radioactive atom distribution image Y, and using K GPUs (GPU 1 to GPU K), the set of radioactive atoms allocated to each GPU is computed in parallel. The allocation relationship is as follows:

[0079]

[0080] in, This represents the set of radioactive atoms allocated to each GPU. In the The value at the individual element, This indicates a round-down operation; it divides the radioactive atom distribution image into K sets of radioactive atoms that have the same total number and spatial distribution.

[0081] d) Following the method in step c), each set of radioactive atoms is assigned to a GPU, and K independent MC simulation tasks are constructed. Each GPU performs an MC simulation task according to the assigned set of radioactive atoms and the attenuation image to perform a preliminary scattering estimation, and obtains a preliminary estimated scattering factor that matches the event.

[0082] 3) Based on the geometric information of PET, the normalized information of detection efficiency, and the attenuation image, the sensitivity image is calculated using the standard process of the EM-type algorithm. If it is necessary to accelerate the sensitivity image calculation process, refer to Chinese Patent CN117648195A (publication date 2024.03.05) to obtain the sensitivity image; steps 1) and 2) are calculated sequentially, and steps 1) and 2) are calculated in parallel and synchronously with step 3).

[0083] 4) Image reconstruction is performed using the EM algorithm. The image reconstruction process includes all physical corrections. Among them, the scattering correction in the physical correction first calculates an activity distribution image with a rough scattering correction based on the initial estimated scattering factor, coincidence event, and sensitivity image. Each iteration calculates an activity distribution image with scattering correction. The scattering factor is then finely estimated using the activity distribution image with scattering correction to obtain the finely estimated scattering factor. The next iteration is then performed. During the iteration, the scattering factor and the activity distribution image become increasingly accurate. The activity distribution image obtained after the iteration converges is the reconstructed image. The fine scattering estimation uses a multi-GPU parallel MC computing framework to perform fine scattering estimation on the scattering factor based on the activity distribution image to obtain the finely estimated scattering factor. The finely estimated scattering factor is then introduced into subsequent iterations to update the activity distribution image, realizing the coordinated iterative update of the scattering factor and the activity distribution image. After 5 iterations, the iterative image converges, and the reconstructed image is obtained.

[0084] The reconstructed images of the NEMA IQ phantom obtained by conventional methods and the reconstructed images obtained by the method of this invention are as follows: Figure 6 As shown, from Figure 6 It can be seen that the quality of the reconstructed image of the present invention is comparable to that of the traditional method; since the reconstruction method of the present invention has one less EM iteration without scattering correction in the overall process than the traditional method, and the present invention optimizes the load distribution of multi-GPU parallel computing, the reconstruction time of the present invention is reduced by about one-third compared with the traditional method.

[0085] Finally, it should be noted that the purpose of disclosing the embodiments is to help further understand the present invention. However, those skilled in the art will understand that various substitutions and modifications are possible without departing from the spirit and scope of the present invention and the appended claims. Therefore, the present invention should not be limited to the content disclosed in the embodiments, and the scope of protection of the present invention is defined by the claims.

Claims

1. A PET scattering correction and image reconstruction method based on time-of-flight and GPU, characterized in that, The method includes the following steps: 1) The positron emission tomography (PET) equipment acquires coincidence events, calculates time-of-flight (TOF) information based on the coincidence events, and calculates the most likely annihilation location based on the geometric information of the PET and the TOF information; the weight of the coincidence event is calculated based on the normalized information of the detection efficiency and the attenuation image; and the TOF histogram image is calculated based on the most likely annihilation location, coincidence event type and weight. 2) The TOF histogram image is used as the initial radiation source distribution for scattering estimation; based on the attenuation image and the initial radiation source distribution, a preliminary scattering estimation of the scattering factor is performed to obtain the preliminary estimated scattering factor that matches the event; in the process of preliminary scattering estimation, atomic-driven multi-GPU parallel computing is used. 3) Based on the geometric information of PET, the normalized information of detection efficiency, and the attenuation image, the sensitivity image is calculated to obtain the sensitivity image; 4) The Expectation-Maximization (EM) algorithm is used for image reconstruction. The image reconstruction process includes all physical corrections. Among them, the scattering correction in the physical correction first calculates the activity distribution image with the initial scattering correction based on the initial estimated scattering factor, coincidence event, and sensitivity image. Then, the scattering factor is finely estimated based on the activity distribution image to obtain the finely estimated scattering factor. The finely estimated scattering factor is introduced into the subsequent iteration to update the activity distribution image. The scattering factor and the activity distribution image are updated in a coordinated iterative manner. The reconstructed image is obtained after the iteration converges.

2. The method according to claim 1, characterized in that, In step 1), the TOF histogram calculation includes the following steps: a) Record matching events: Each coincidence event includes: the crystal position of the fast photon. Fast Photon Timestamp Crystal position of slow photons Slow photon timestamp and matching event type The event types include instantaneous events and delayed events. For instantaneous events... =0, for delayed events =1; b) Calculate TOF information: Read the timestamps of two single-photon events from the matching events, and calculate the Time-of-Flight (TOF) information based on the timestamps of the two single-photon events. for: c) Calculate the most likely annihilation location: Based on the TOF information and the geometric information of the crystal positions of the two single-photon events (i.e., PET), the most likely annihilation location is calculated as follows: Where c represents the speed of light. The most likely annihilation location; d) Calculate the weights: Based on the normalized information of PET detection efficiency, the normalized correction factor corresponding to the matching event is calculated as follows: Perform orthographic projection calculations on the attenuation image to determine the attenuation correction factor corresponding to the event. Then calculate the weight corresponding to the matching event. for: e) Calculate voxel values: Given a TOF histogram image X, where all voxels are initially set to 0, update the voxel values ​​of the TOF histogram image X using all matching events; if the most probable annihilation location... The first point that falls on the TOF histogram image X In individual elements, then according to the event... and weight Update # Individual physiology value : Where j=1,..,J, J is the number of voxels in the TOF histogram image X until all coincidence events have been calculated, thus obtaining the TOF histogram image.

3. The method according to claim 1, characterized in that, In step 2), the Monte Carlo simulation (MC) method or the single scattering event simulation (SSS) method is used to perform a preliminary scattering estimate of the scattering factor.

4. The method according to claim 3, characterized in that, The atomic-driven multi-GPU parallel MC computing framework includes the following steps: a) Convert the TOF histogram image X into a radioactive atom distribution image Y. The conversion relationship is as follows: in, and The first two images represent the TOF histogram image X and the radioactive atom distribution image Y, respectively. Individual phenotypic value The total number of radioactive atoms to be simulated; b) Let the PET acquisition time simulated by the MC method be... The total number of radioactive atomic decays simulated by the MC method is M: in, The half-life of a radioactive nuclide; c) Employing an atomic-level load balancing strategy, the set of radioactive atoms allocated to each GPU is calculated through multi-GPU parallel computing with a total of K GPUs, by traversing all voxels of the radioactive atom distribution image Y. The allocation relationship is as follows: in, This represents the set of radioactive atoms allocated to each GPU. In the The value at the individual element, This indicates a round-down operation; d) Each GPU performs a preliminary scattering estimate based on the attenuation image according to the assigned set of radioactive atoms, and obtains a preliminary estimated scattering factor that matches the event.

5. The method according to claim 1, characterized in that, In step 4), each iteration calculates an activity distribution image with scattering correction. The scattering factor is then estimated using the activity distribution image with scattering correction to obtain the fine scattering factor, and the next iteration is performed. The activity distribution image obtained after the iteration converges is the reconstructed image.

6. The method according to claim 5, characterized in that, Detailed scattering estimation employs a multi-GPU parallel MC computation framework.

Images