GPU-accelerated three-dimensional reconstruction method of the Compton camera and its application in radionuclide imaging

Ren-Yao Wu ·Chang-Ran Geng ·Feng Tian ·Zhi-Yang Yao ·Chun-Hui Gong ·Hao-Nan Han ·Jian-Feng Xu,4 ·Yong-Shun Xiao ·Xiao-Bin Tang

Abstract A novel and fast three-dimensional reconstruction method for a Compton camera and its performance in radionuclide imaging is proposed and analyzed in this study.The conical surface sampling back-projection method with scattering angle correction (CSS-BP-SC) can quickly perform the back-projection process of the Compton cone and can be used to precompute the list-mode maximum likelihood expectation maximization (LM-MLEM).A dedicated parallel architecture was designed for the graphics processing unit acceleration of the back-projection and iteration stage of the CSS-BP-SC-based LM-MLEM.The imaging results of the two-point source Monte Carlo (MC) simulation demonstrate that by analyzing the full width at half maximum along the three coordinate axes, the CSS-BP-SC-based LM-MLEM can obtain imaging results comparable to those of the traditional reconstruction algorithm, that is, the simple back-projection-based LM-MLEM.The imaging results of the mouse phantom MC simulation and experiment demonstrate that the reconstruction results obtained by the proposed method sufficiently coincide with the set radioactivity distribution, and the speed increased by more than 664 times compared to the traditional reconstruction algorithm in the mouse phantom experiment.The proposed method will further advance the imaging applications of Compton cameras.

Keywords Compton camera ·Three-dimensional reconstruction ·Radionuclide imaging ·GPU

1 Introduction

Radionuclide imaging allows the noninvasive visualization of the distribution of radionuclides in the region of interest and is an important step in the diagnosis and treatment of several diseases, including cancer [1, 2, 3].Single-photon emission computed tomography (SPECT) and positron emission computed tomography (PET) are currently the two widely used radionuclide imaging techniques.However, the energy ranges of gamma rays that can be detected are limited to 511 keV (PET) or generally below 300 keV(SPECT) [4].Considering the aforementioned, the Compton camera (CC) based on Compton scattering kinematics has been proposed for medical imaging owing to its wide detection energy range, multi-radionuclide detection capability,and high detection efficiency [5, 6, 7, 8, 9].Studies indicate that a CC presents significant potential to be used as nextgeneration detection equipment [10, 11, 12].

Todd et al.first introduced the CC concept to nuclear medical imaging in 1974 [13].Although capturing low-energy gamma rays is difficult for a CC owing to the energy dependence of the probability of Compton scattering occurring, several groups have demonstrated CC imaging results for99mTc (141 keV), a radionuclide commonly used in traditional SPECT [14, 15].In recent years, the application of CCs in in vivo radionuclide imaging and dose monitoring of prompt gamma rays has been extensively studied[10, 11, 12, 13, 14, 15, 16, 17, 18]; Takashi et al.performed a human clinical trial of in vivo radionuclide imaging based on a CC in 2020 [12].The CC imaging results are related to a series of factors, such as the detector material, energy resolution, and angular resolution.Generally, semiconductor-based CCs are widely developed owing to their excellent detection performance [19].Figure 1 illustrates the principle of the CC system.Based on the recorded scattering and absorption information of Compton event pairs, a CC can typically locate the radioactive source on a conical surface with the scattering location as its vertex and the Compton scattering angle as its apex angle [20, 21].Multiple Compton event pairs produce multiple conical surfaces, and the intersection of these conical surfaces is the actual location of the source [22].

The image reconstruction of CC relies on the back-projection of conical surfaces into an imaging space, which is difficult to represent mathematically using an algorithm.Several studies have indicated that image reconstruction from CC detection data is challenging [23, 24, 25].So far,the image reconstruction methods of a CC include direct analytical and statistical iterative methods [26, 27, 28, 29,30].Simple back-projection (SBP) is a direct conventional CC image reconstruction method in which the voxelized imaging space is traversed, and all the voxels that intersect the conical surfaces generated by the Compton events within the imaging space are found.Therefore, the SBP result is a collection of sets of voxels associated with each Compton conical surface in a voxelized imaging space [31].However,the reconstruction results of SBP generally suffer from poor image resolution, a low signal-to-noise ratio, and apparent artifacts.Therefore, the maximum likelihood expectation maximization (MLEM) iterative algorithm is recommended for high-resolution reconstructions [32].MLEM allows the approximation of images of radioactive sources by iteratively sharpening a coarse approximation, where list-mode MLEM (LM-MLEM) is a currently used iterative form of the optimized MLEM [33, 34].Studies have also proposed using the Markov-chain-based stochastic origin ensemble(SOE) algorithm for CC image reconstruction in nuclear medical imaging [35, 36, 37].For medical imaging, the enhanced speed of image reconstruction translates directly into clinical workflow benefits [38].The rapid imaging of radionuclides facilitates timely feedback of clinical information, such as supporting intraoperative tumor localization and image-guided interventions [39, 40].However, owing to the complex iterative calculation of conical back-projection and three-dimensional (3D) reconstruction, the reconstruction speed of the aforementioned CC image reconstruction methods makes it difficult to meet clinical needs.

Graphics processing units (GPUs) have been previously used in studies to accelerate the reconstruction process of the CC to overcome the challenge of reconstruction speed.Zheng et al.performed GPU acceleration on the SOE algorithm, where each Compton event detected by the CC corresponded to each thread on the GPU and achieved a 25 × increase in speed compared to the serial program[41].However, the parallel SOE algorithm can be further improved in terms of convergence and 3D reconstruction.Nauyen et al.performed parallel acceleration on the ordered subset expectation maximization algorithm based on the ray tracing method and achieved a 16 × increase in speed compared to the serial program [42].The key to this approachis storing the rays forming a conic surface, which is computationally complex and requires large memory.Yao et al.proposed a rapid subset-driven origin ensemble with a resolution recovery (SD-OE-RR) algorithm that can correct the image reconstruction bias caused by the uncertainty of the detection data of the CC [43]; however, this did not accelerate the iterative part of the algorithm.

Fig.1 Block diagram and principle of the CC system

In this study, we propose a novel GPU-accelerated CC image reconstruction method.A conical surface sampling method based on a coordinate system transformation was applied to reduce the computational complexity of the Compton cone-based back-projection process.A Compton scattering angle correction method was introduced in the conical surface sampling process.The reconstructed results of the Compton cone-based back-projection served as the pre-calculation of the LM-MLEM to perform the iterative process.In addition, through a dedicated parallel architecture designed for the GPU acceleration of the backprojection and LM-MLEM iteration stage, our method significantly improves the 3D reconstruction speed of the CC.

2 Three-dimensional reconstruction method

2.1 Back-projection of the Compton cone

2.1.1 Conical surface sampling method based on the transformation of the coordinate system

In this study, a 3D position-sensitive CdZnTe CC (3DCZT CC) was used for the CC image reconstruction.As shown in Fig.2a, a circular 3D-CZT CC array collects the Compton event data from different angles.Through the transformation of the coordinate system, a flexible expression of the conical surface in a 3D space can be achieved.Fast image reconstruction is achieved by uniform spatial sampling on the conical surfaces.

As shown in Fig.2b,αandβare the angles formed by the projection of -→uin theXZplane and theZ-axis, and the-→uitself.The laboratory coordinate system can be transformed using the rotation matrixMas follows:

The transformation results are shown in Fig.2c.The quadratic equation of a conical surface in this new coordinate system can be expressed by Eq.(3), indicating a standardized cone in a mathematical expression:

The transformed coordinate system is called the Compton event coordinate system.The vertex coordinates of the conical surface in this coordinate system can be described as follows:

Fig.2 (Color online) Schematic of the Compton event reconstruction.a Schematic diagram of multi-angle detection.b Compton back-projection cone in the laboratory coordinate system.c Compton back-projection cone in the Compton event coordinate system

After the rotation of the coordinate system, the sampling space for the process of conical surface sampling in the Compton event coordinate system should be a circumscribed sphere of the imaging space, as shown in Fig.2c.

Meanwhile, (x�,y�,z�) are the coordinates on the conical surface.In the conical surface sampling process,z′can be determined in the direction of the axis of the standardized cone, and its parameter can be determined by uniform sampling as follows:

In the conical surface sampling process, each voxel at the same location in the imaging space is only recorded once for each Compton event.After (x�,y�,z�) is obtained, an inverse transformation can be applied to obtain the coordinates of(x,y,z) in the laboratory coordinate system.

2.1.2 Correction of the Compton scattering angle

The 3D-CZT CC prototype used in this study is fabricated by the Kromek group [44, 45].The size of the CZT crystal was 22 mm × 22 mm × 15 mm, which was divided into 11 × 11 pixelated anodes with one planar cathode; each pixel area was 2 mm × 2 mm.The spatial resolution over the entire 3D-CZT volume was 2.0 mm × 2.0 mm × 0.34 mm.The energy resolutionδE(full width at half maximum (FWHM)) of the 3D-CZT CC at different energy levels was fitted through the experimental results, as shown in Fig.3a.

The Compton scattering angleθis calculated using Eq.(7):

wheremec2is the electron rest mass,E0is the incident gamma-ray energy,E1is the measured deposition energy of the recoil electrons by Compton scattering, andE2is the measured deposition energy of the scattered photons.In radionuclide imaging, the energy of the characteristic gamma rays produced by the radionuclide is generally known; therefore, the Compton scattering angle error caused by the uncertainty of the measured energy should be considered.In this study, the influence of the energy resolution and Doppler effect on the energy detected by 3D-CZT was modeled as an angular Gaussian blur of the Compton cone angleθ[46, 47].The Gaussian blurring method for the cone angle can be found in Eq.(8) as follows:

Fig.3 (Color online) Correction of conical angle.a Energy resolution of the 3D-CZT CC.b Correction of the Compton scattering angle θ

whereσθrepresents the sigma of angular Gaussian blurring.For the conical surface sampling process described in the previous section, a newθ′is first generated according to Eq.(8) and subsequently introduced into Eqs.(3)-(6); this process is continuously repeated.Figure 3b presents the geometric relationship of Eq.(8), withσθindividually calculated for each Compton event.This process guarantees a Gaussian distribution of the sampling points around the cone angleθ.

The aforementioned image reconstruction method is referred to as the conical surface sampling back-projection with scattering angle correction (CSS-BP-SC) algorithm,which includes the conical surface sampling back-projection (CSS-BP) and Compton scattering angle correction.In this study, when a voxel is sampled, a weight of 1 is automatically increased to that of the voxel.

2.2 Implementation of the MLEM iteration

The LM-MLEM iterative algorithm is used to image the distribution of the radioactive sources.The iterative process of the LM-MLEM algorithm is expressed by the following equation:

where,jis the index of voxels in the imaging space,iis the index of detected Compton events,fl

jis the intensity of thejth voxel in the imaging space afterliterations,Sjis the relative probability of detecting the Compton events originating from thejth voxel, and the system matrixtijcorresponds to the probability of a photon emitted by thejth voxel to be detected as theith Compton event.

Whenlis equal to zero, the distribution offl jis the result of the Compton conical back-projection, that is, the precomputation of the LM-MLEM iteration.Moreover,Sjwas set to be uniform throughout the imaging space in this study [48].The analytical calculation of the system matrix has been discussed by Maxim et al.;tijis finally positively related to the differential cross section of the Compton scattering for theith Compton event [34].

The forward projection of the LM-MLEM iteration contributes to the likelihood estimate of each Compton event:

The forward projection of the LM-MLEM iteration is the summation of the system matrix element of each Compton eventimultiplied by the image distribution.Thetikdepends only on indexiand is nonzero on the set of voxels that contribute to theith Compton event; thus, the following equations can be obtained:

wheretirepresents the factor positively related to the differential cross section of theith Compton event, andFlirepresents the sum of the intensities of the set of voxels that contribute to theith Compton event.Based on the forward projection, the back-projection process of the LM-MLEM iteration can be expressed as follows:

The back-projection process of the LM-MLEM iteration is the summation of the ratio of the system matrix element of each Compton evention voxeljto its forward projection,also called the back-projection factor of thejth voxel.In the calculation of the back-projection factor,tijrepresents the system matrix element for Compton eventi, whose value is only related to indexi; therefore, the factor can be eliminated along withtiin the denominator.Subsequently, the backprojection factor can be calculated as follows:

The proposed 3D reconstruction method is referred to as the CSS-BP-SC-based LM-MLEM iterative reconstruction algorithm, which uses CSS-BP-SC as the precomputation of the iteration.In this study, the imaging space was divided into 128 × 128 × 128 voxels.The side length of each cubic voxel was approximately 1.09 mm, and the conical surface sampling number was set at 2.4 × 105.

2.3 Parallel computing architecture for reconstruction

This section presents the implementation of the proposed CSS-BP-SC-based LM-MLEM iterative reconstruction algorithm in a parallel architecture.Compute unified device architecture (CUDA) is used in the parallel computing of the CSS-BP-SC-based LM-MLEM reconstruction algorithm.Parallel computation can be divided into the following three parts: CSS-BP-SC, LM-MLEM iteration parameters,and LM-MLEM iterative process.Parallel computing is supported by GPU hardware.

1.Parallel computation of CSS-BP-SC:

For each detected event, the CSS-BP-SC process is performed by a thread on the GPU, which generates the corresponding probability density estimation (PDE).Simultaneously, the total density matrix is obtained by summing the PDE of each thread.The total density matrix describes the distribution off0j.The total number of effective samplings ofLifor theith detected event is recorded.

2.Parallel computation of the LM-MLEM iteration parameters:

The collection of all the voxels intersecting the conical surfaces produced by each Compton event in the imaging space is called a voxel set for the Compton event.The key to Eqs.(12) and (14) is establishing the mapping relationship between each Compton event and voxel in its associated voxel set.The index of each voxel in each voxel set is referred to as the LM-MLEM iteration parameter.Considering that the memory address available for GPU operations must be continuous, a parameter matrix having a lengthL(L=L1+L2+…+LI) is constructed in the video memory to store the LM-MLEM iteration parameters, whereIis the total number of detected events.Therefore, each Compton event corresponds to part of the parameter matrix.

3.Parallel computation of the LM-MLEM iterative process:

A simple description of the parallel architecture is presented in Fig.4.The parallel iterative reconstruction algorithm was programmed using the hybrid C++ and CUDA C languages.In this study, GPU hardware parallel acceleration was implemented using an NVIDIA Titan V with CUDA 11.4.In contrast, the reconstruction algorithms without GPU acceleration were executed with one thread on a single core on the Linux platform using an Intel Xeon processor E5-2699 v4.

3 Simulation and experimental settings

We performed both the Monte Carlo (MC) simulations and experiments to evaluate the performance of the proposed algorithm.

3.1 Description of the Monte Carlo simulations

The MC toolkit Geant4 (version 10.05) was used for the simulations, and the G4EmPenelopePhysics module was added to the built-in physics list QGSP_BIC of Geant4 for precise photon transportation.As shown in Fig.2a, a multiangle CC array composed of 3D-CZTs was constructed.Two 511 keV mono-energy point sources were simulated,as shown in Fig.5a.The activity of the two-point source was set to be identical.Furthermore, the practicability of the proposed algorithm for radionuclide imaging applications was evaluated using a mouse phantom for simulation.The geometry of the simulation is shown in Fig.5b.The mouse phantom was referenced from a normal 28 g male mouse and generated by co-registered CT and cryosection images [49].In the simulation, the phantom was divided into 380 × 992 × 208 voxels, and each voxel was 1 × 10-3mm3in size.The mouse phantom contained the following four organs: the brain, heart, kidneys, and bladder.The brain and bladder were analyzed in this study.The material of the torso of the mouse phantom was butanediol dimethacrylate(C12H18O4,ρ: 1.3 g/cm3), while that of the organs of the mouse was water.Gamma rays of 511 keV are emitted from the brain and bladder in the mouse phantom simulation and have the same intensity.

Fig.4 Architecture of parallel computing

Fig.5 (Color online) Schematic of MC simulations.a Side view of the two-point source simulation.b Top view of the mouse phantom simulation

Fig.6 (Color online) Mouse phantom experiment.a Compton camera imaging system.b Mouse phantom experiment.c Electric turntable.d Illustration of the geometry and radiation source settings of the experiment

The MC simulation results can provide an accurate energy deposition, which is impossible in practical situations.Therefore, in the post-processing of the Compton data,Gaussian broadening equivalent to the actual energy resolution was applied to the recorded energies of each Compton event pair to adapt to the actual detection situation; the actual energy resolution was derived from the fitted values shown in Fig.3a.In both the simulations and experiments in this study, we selected an energy window of ± 5 keV to screen for valid Compton events.As shown in Fig.5, nine two-point source and mouse phantom simulations were performed with the detector array rotated (10° per rotation)to detect the 36 angles.The Geant4 codes were executed by a 64-bit Linux computer using an Intel Xeon processor E5-2699 v4.

3.2 Description of the mouse phantom experiment

Based on the simulations introduced in the previous section,we conducted a mouse phantom experiment with the same geometry and radiation source setup, as shown in Fig.6.The 3D-printed mouse phantom used in the experiment was identical to that used in the simulation.In the experiment,the mouse phantom was driven by an electric turntable for multi-angle detection (Qingdao Shenji Intelligent Technology Co., Ltd.; better than 0.1° accuracy).A total of 36 angles were detected, and the detection time for each angle was 45 s.The radiation source was a 511 keV gamma-ray emitter[18F] sodium fluoride ([18F]NaF), which was produced by JYAMS PET Research and Development Limited.The [18F]NaF solution was sealed in capsules and placed at the brain and bladder sites of the mouse phantom with a radioactivity of ～ 25.1 and ～ 24.7 μCi, respectively.Referring to the study by Shy et al.for the ordering rule of the Compton event pairs in the experiment, with an incident gamma-ray energy above 350 keV, the energy deposition of scattering events is generally considered greater than that of absorption events [50].

4 Results

4.1 Reconstruction results of the two-point source simulation

The reconstruction results for the two-point source were first evaluated using the FWHM along the three coordinate axes.Figure 7 presents the reconstruction results of using SBP,CSS-BP, and CSS-BP-SC.The FWHM of the SBP, CSS-BP,and CSS-BP-SC results is 15.53, 15.19, and 15.02 mm along thex-axis 22.16, 22.14, and 22.46 mm along they-axis, and 22.42, 22.47, and 22.46 mm along thez-axis, respectively;the reconstruction times of these three methods are 16,046.5,1117.2, and 1223.3 s, respectively.Figure 8 presents the LMMLEM iterative reconstruction results using the SBP, CSS-BP,and CSS-BP-SC methods for precomputation.The FWHM of the SBP-based, CSS-BP-based, and CSS-BP-SC-based LMMLEM is 7.22, 7.24, and 7.39 mm along thex-axis, 9.17, 9.24,and 9.42 mm along they-axis, and 9.01, 8.93, and 9.19 mm along thez-axis, respectively; the reconstruction times of these three methods are 16,339.7, 1357.9, and 1559.5 s, respectively.The reconstruction results shown in Figs.7 and 8 are the slices containing the radioactive source in the corresponding 3D reconstruction results.

The calculation of the contrast (Vcon) is expressed by Eq.(16),

Fig.7 (Color online) Reconstruction results of SBP, CSS-BP, and CSS-BP-SC.Reconstruction images by a SBP, b CSS-BP, and c CSS-BP-SC.Probability density distribution along the d x-axis, e y-axis, and f z-axis

where MeanROIis the average pixel intensity of the region of interest (ROI), and the MeanBackgroundis the average pixel intensity in the background.The ROI is determined by using the axial FWHM value of the reconstruction results.For the iterative reconstruction results shown in Fig.8, the calculatedVconvalues of the SBP-based, CSS-BP-based, and CSS-BP-SC-based LM-MLEM are 282.08, 267.94, and 292.73, respectively.A total of 30,676 Compton events were used in the reconstruction.The FWHM along the three axes and the reconstruction times for the aforementioned methods are listed in Table 1.

4.2 Reconstruction results of the mouse phantom simulation

Figure 9 presents the iterative reconstruction results of a mouse phantom in the simulation, which were obtained using different methods: SBP-based, CSS-BP-SC-based, and GPU-accelerated CSS-BP-SC-based LM-MLEM methods.A total of 34,539 Compton events were used in the reconstruction.The reconstructed images of these three methods matched the real distribution.Mouse phantom reconstruction results were quantified using the activity recovery coefficient (ARC).For each radioactive organ, ARC is defined as the ratio between the activity reconstructed in the organ position and its true activity as follows [51]:

whereAiandViindicate the activity reconstructed in radioactive organiand its volume, respectively, andATis the total reconstructed activity in the entire image, either inside or outside the radioactive organ volume.The brain and bladder are radiopharmaceutical-enriched organs.For the brain,the ARCs of the SBP-based, CSS-BP-SC-based, and GPUaccelerated CSS-BP-SC-based LM-MLEM reconstruction results were 0.243, 0.240, and 0.240, respectively.For the bladder, the ARCs of these three reconstruction methods were 0.294, 0.295, and 0.294, respectively.Table 2 lists the time consumption of the mouse phantom reconstruction process based on these three methods.

4.3 Reconstruction results of the mouse phantom experiment

Fig.8 (Color online) Reconstruction results of the iterative reconstruction methods.Reconstruction images of the a SBP-based, b CSS-BPbased, and c CSS-BP-SC-based LM-MLEM methods.Probability density distribution along the d x-axis, e y-axis, and f z-axis

Figure 10 presents the reconstruction results of the mouse phantom used in the experiment.A total of 20,271 Compton events were used for reconstruction.The timeconsumptions of the SBP-based, CSS-BP-SC-based, and GPU-accelerated CSS-BP-SC-based LM-MLEM methods are 10,750.7, 1102.8, and 16.2 s, respectively.The radiopharmaceutical is assumed to be uniformly distributed throughout the target organs by default; therefore,the radioactivity areas in the mouse phantom experiment are also considered to be associated with the real positions of the organs.For the brain, the ARCs of the SBPbased, CSS-BP-SC-based, and GPU-accelerated CSS-BPSC-based LM-MLEM reconstruction results were 0.047,0.049, and 0.047, respectively.For the bladder, the ARCs of these three reconstruction methods were 0.037, 0.035,and 0.034, respectively.The GPU-accelerated reconstruction method was more than 664 times faster than the traditional reconstruction algorithm, that is, the SBP-based LM-MLEM.Although the experimental results of the radionuclide imaging in the mouse phantom were inferior to those of the simulation, the radiation sources from the brain and bladder were clearly identified.Table 3 lists the time consumption of the mouse phantom reconstruction process based on these three methods.

Table 1 Reconstruction performance of the different methods for the reconstruction of the two-point source

5 Discussion

A novel GPU-accelerated CSS-BP-SC-based LM-MLEM iterative reconstruction method is proposed, and the performance of the algorithm applied to radionuclide imaging is studied using 3D-CZT CC.

Considering the 3D reconstruction method, Fig.7 demonstrates that the SBP, CSS-BP, and CSS-BP-SC methods have comparable imaging results, while the reconstruction speeds of the CSS-BP and CSS-BP-SC methods can reach 14 times and 13 times that of the SBP, respectively.Owing to the correction of the scattering angle, Fig.8 demonstrates that the CSS-BP-SC-based LM-MLEM improves the contrastVconby 9.3% compared to that of the CSSBP-based LM-MLEM.Owing to the acceleration of the back-projection stage, the reconstruction speed of the CSS-BP-SC-based LM-MLEM can reach 11 times that of the SBP-based LM-MLEM in the reconstruction of the two-point source.The high-speed reconstruction of the CSS-BP-SC in the back-projection stage benefits from the fact that the sampling dimension of the spatial sampling method is independent of the reconstructed spatial parameters.In particular, when the reconstruction space is finely divided, the projection process of SBP must traverse more voxels unrelated to the back-projection Compton cone,which is not required for CSS-BP-SC.

Fig.9 (Color online) Reconstruction results of the SBP-based, CSSBP-SC-based, and GPU-accelerated CSS-BP-SC-based LM-MLEM methods of the mouse phantom in simulation.Each figure indicates that the slices are 2.18 mm apart from one another in the direction perpendicular to the mouse coronal plane

Table 2 Time consumption of the different reconstruction methods with the simulation data

Fig.10 (Color online) Reconstruction results of the SBP-based, CSSBP-SC-based, and GPU-accelerated CSS-BP-SC-based LM-MLEM methods of the mouse phantom in the experiment.Each figure indicates that the slices are 2.18 mm apart from one another in the direction perpendicular to the mouse coronal plane

Table 3 Time consumption of the different reconstruction methods with the experimental data

In this study, we designed a specialized parallel architecture for CSS-BP-SC-based LM-MLEM, as shown in Fig.4.In the mouse phantom simulation, Fig.9 and Table 2 demonstrate the reconstruction acceleration effect of the proposed 3D iterative reconstruction method under GPU hardware parallel acceleration.Compared to the serial CSS-BP-SC-based LM-MLEM, the parallel-executed CSSBP-SC-based LM-MLEM achieves an increase in speed of greater than 67 times in 3D reconstruction.The slight difference in the reconstruction results between the accelerated and unaccelerated programs arises from the difference in the random-number seed used in the CSS-BP-SC process.Meanwhile, the GPU-accelerated CSS-BP-SC-based LM-MLEM iterative method demonstrates an increase in speed of greater than 722 times compared to the SBP-based LM-MLEM iterative method.The video memory size is the limiting factor for the GPU acceleration ratio; in multi-GPU parallelism, the reconstruction speed will also further increase linearly.The GPU-accelerated CSS-BP-SC-based LM-MLEM iterative method demonstrated similar evaluated ARC values compared to the SBP-based LM-MLEM iterative method, as well as similar to the values in a previous study [51].

In the mouse phantom experiment, the mouse phantom was detected at 36 angles, each with a detection time of 45 s.Therefore, in a real ring detector, radionuclide detection requires only a few minutes.Regarding the reconstruction results, the proposed GPU-accelerated method demonstrates ARC evaluation values comparable to those of the traditional reconstruction algorithm.We found that the experimental results for the mouse phantom were inferior to the simulation results.This may be due to the electronic noise or effects of the adjacent pixels.These limitations are related to the properties of CC.In addition, the Compton event ordering rule used in the experiment is one of the reasons for the deterioration of the reconstruction results [50].When the incident ray has an energy of 511 keV, the probability of the energy of the scattering in the Compton event pair being greater than the energy of the absorption is less than 60% [52].Therefore,the sequence of events in a considerable portion of the experimental results would be incorrect and thus become artifacts in the reconstructed images.The radionuclide imaging results of the mouse phantom experiment are comparable to other CC-based radionuclide imaging studies [8, 9, 10, 11, 12, 23], providing a basis for the application of our reconstruction method to the 3D radionuclide imaging of 3D-CZT CC.Although we used 3D-CZT as the detector, virtually any CC currently known is suitable for our proposed 3D reconstruction method.The proposed CC imaging method is suitable for the range detection of proton therapy or monitoring of the boron concentration in boron neutron capture therapy (BNCT), in addition to radionuclide imaging [53, 54].In particular, rapid imaging of the particle range information in proton therapy and boron concentration distribution information in BNCT therapy.

However, certain limitations should be considered in future studies.In follow-up studies, an accurate and rigorous system and sensitive matrices can be applied to the LM-MLEM iteration.The setup of the radioactive source in the current phantom experiment was relatively simple,and a phantom experiment composed of several hot spots surrounded by a warm background can be performed in the future.Imaging results from multi-angle detection should be compared with existing nuclear medicine imaging devices,such as SPECT or PET; animal studies are also necessary for the preclinical validation of this reconstruction method.Only a 511 keV gamma source was analyzed in this study;however, simultaneous imaging experiments of radionuclides with multiple energy levels can extend the application of the reconstruction method in nuclear medicine imaging.

6 Conclusion

In this study, we propose the use of a GPU-accelerated CSSBP-SC-based LM-MLEM iterative method for 3D imaging with a CC.The proposed GPU-accelerated method dramatically improves the computational performance of the CC image reconstruction while demonstrating imaging results comparable to those of traditional reconstruction algorithms.The effectiveness and accuracy of the method were validated through MC simulations and experiments using a mouse phantom.These results demonstrate the promising applications of the proposed methodology.

Author contributionsAll authors contributed to the study conception and design.Material preparation, data collection, and analysis were performed by Ren-Yao Wu, Chang-Ran Geng, Feng Tian, Zhi-Yang Yao, Chun-Hui Gong, Hao-Nan Han, Jian-Feng Xu, Yong-Shun Xiao, Xiao-Bin Tang.The first draft of the manuscript was written by Ren-Yao Wu, and all authors commented on previous versions of the manuscript.All authors read and approved the final manuscript.

Data availabilityThe data that support the findings of this study are openly available in Science Data Bank at https:// doi.org/ 10.57760/ scien cedb.07763 and https:// cstr.cn/ 31253.11.scien cedb.07763.