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

GPU-accelerated three-dimensional reconstruction method of the Compton camera and its application in radionuclide imaging 1 first-author 1 corresp * gengchr@nuaa.edu.cn * gengchr@nuaa.edu.cn gengchr@nuaa.edu.cn https://orcid.org/0000-0002-6332-6193 1 2 https://orcid.org/0000-0001-7344-4460 3 1 1 4 1 1 corresp tangxiaobin@nuaa.edu.cn tangxiaobin@nuaa.edu.cn tangxiaobin@nuaa.edu.cn https://orcid.org/0000-0003-3308-0468 Department of Nuclear Science and Technology, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, China Department of Nuclear Science and Technology, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, China Department of Engineering Physics, Tsinghua University, Beijing 100084, China Department of Engineering Physics, Tsinghua University, Beijing 100084, China School of Environmental and Biological Engineering, Nanjing University of Science and Technology, Nanjing 210094, China School of Environmental and Biological Engineering, Nanjing University of Science and Technology, Nanjing 210094, China JYAMS PET Research and Development Limited, Nanjing 211100, China JYAMS PET Research and Development Limited, Nanjing 211100, China 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. A 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. Compton camera Three-dimensional reconstruction Radionuclide imaging GPU Introduction 1 Radionuclide imaging allows the non-invasive 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 1 , 2 2 , 3 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 4 ]. Considering the aforementioned, a 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 5 , 6 6 , 7 7 , 8 8 , 9 9 ]. Studies indicate that a CC presents significant potential to be used as next-generation detection equipment [ 10 10 , 11 11 , 12 12 ]. Todd et al. first introduced the CC concept to nuclear medical imaging in 1974 [ 13 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 for 99m Tc (141 keV), a radionuclide commonly used in traditional SPECT [ 14 14 , 15 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 10 , 11 11 , 12 12 , 13 13 , 14 14 , 15 15 , 16 16 , 17 17 , 18 18 ]; Takashi et al. performed a human clinical trial of in vivo radionuclide imaging based on a CC in 2020 [ 12 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 19 ]. Fig. 1 Fig. 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 20 , 21 21 ]. Multiple Compton event pairs produce multiple conical surfaces, and the intersection of these conical surfaces is the actual location of the source [ 22 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 23 , 24 24 , 25 25 ]. Thus far, the image reconstruction methods of a CC include direct analytical and statistical iterative methods [ 26 26 , 27 27 , 28 28 , 29 29 , 30 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 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 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 33 , 34 34 ]. Studies have also proposed using the Markov-chain-based stochastic origin ensemble (SOE) algorithm for CC image reconstruction in nuclear medical imaging [ 35 35 , 36 36 , 37 37 ]. For medical imaging, the enhanced speed of image reconstruction translates directly into clinical workflow benefits [ 38 38 ]. The rapid imaging of radionuclides facilitates timely feedback of clinical information, such as supporting intraoper ative tumor localization and image-guided interventions [ 39 39 , 40 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 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 42 ]. The key to this approach is 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 43 ]; however, this did not accelerate the iterative part of the algorithm. 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 back projection and LM-MLEM iteration stage, our method significantly improves the 3D reconstruction speed of the CC. Three-dimensional reconstruction method 1 Back-projection of the Compton cone 2 Conical surface sampling method based on the transformation of the coordinate system 3 In this study, a 3D position-sensitive CdZnTe CC (3D-CZT CC) was used for the CC image reconstruction. As shown in Fig. 2a 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. Equation (1) is the quadratic equation describing the Compton back-projection conical surface [ 23 23 ], where <math id="M1"><mi>θ</mi></math> is the Compton scattering angle, <math id="M2"><mover accent="true"><mi>u</mi><mo>→</mo></mover></math> is the unit vector in the direction of the Compton axis, and <math id="M3"><mrow><mrow><mo>(</mo><mrow><mi>x</mi><mo>−</mo><msub><mi>x</mi><mi>s</mi></msub><mo>,</mo><mtext> </mtext><mi>y</mi><mo>−</mo><msub><mi>y</mi><mi>s</mi></msub><mo>,</mo><mtext> </mtext><mi>z</mi><mo>−</mo><msub><mi>z</mi><mi>s</mi></msub></mrow><mo>)</mo></mrow></mrow></math> represent the vectors connecting the points ( x , y , z ) on the conical surface and the Compton scattering points ( x s , y s , z s ), <math display="block" id="M4"><mtable columnalign="left"><mtr><mtd><msup><mrow><mo>(</mo><mrow><mover accent="true"><mi>u</mi><mo>→</mo></mover><mrow><mo>(</mo><mrow><mi>x</mi><mo>−</mo><msub><mi>x</mi><mi>s</mi></msub><mo>,</mo><mi>y</mi><mo>−</mo><msub><mi>y</mi><mi>s</mi></msub><mo>,</mo><mi>z</mi><mo>−</mo><msub><mi>z</mi><mi>s</mi></msub></mrow><mo>)</mo></mrow></mrow><mo>)</mo></mrow><mn>2</mn></msup><mo>−</mo><mrow><mo>(</mo><mrow><mrow><mo>(</mo><mrow><mi>x</mi><mo>−</mo><msub><mi>x</mi><mi>s</mi></msub><mo>,</mo><mi>y</mi><mo>−</mo><msub><mi>y</mi><mi>s</mi></msub><mo>,</mo><mi>z</mi><mo>−</mo><msub><mi>z</mi><mi>s</mi></msub></mrow><mo>)</mo></mrow></mrow></mrow></mtd></mtr><mtr><mtd><mo>⋅</mo><mrow><mrow><mrow><mo>(</mo><mrow><mi>x</mi><mo>−</mo><msub><mi>x</mi><mi>s</mi></msub><mo>,</mo><mi>y</mi><mo>−</mo><msub><mi>y</mi><mi>s</mi></msub><mo>,</mo><mi>z</mi><mo>−</mo><msub><mi>z</mi><mi>s</mi></msub></mrow><mo>)</mo></mrow></mrow><mo>)</mo></mrow><msup><mtext>cos</mtext><mn>2</mn></msup><mi>θ</mi><mo>=</mo><mn>0.</mn></mtd></mtr></mtable></math> As shown in Fig. 2b Fig. 2b , <math id="M5"><mi>α</mi></math> and <math id="M6"><mi>β</mi></math> are the angles formed by the projection of <math id="M7"><mrow><mo>−</mo><mover accent="true"><mi>u</mi><mo>→</mo></mover></mrow></math> in the XZ plane and the Z -axis, and the <math id="M8"><mrow><mo>−</mo><mover accent="true"><mi>u</mi><mo>→</mo></mover></mrow></math> itself. The laboratory coordinate system can be transformed using the rotation matrix <math id="M9"><mi>M</mi></math> as follows: <math display="block" id="M10"><mrow><mi>M</mi><mo>=</mo><mrow><mo>[</mo><mrow><mtable><mtr><mtd><mrow><mtext>cos</mtext><mi>α</mi></mrow></mtd><mtd><mn>0</mn></mtd><mtd><mrow><mo>−</mo><mtext>sin</mtext><mi>α</mi></mrow></mtd></mtr><mtr><mtd><mrow><mo>−</mo><mtext>sin</mtext><mi>α</mi><mtext>sin</mtext><mi>β</mi></mrow></mtd><mtd><mrow><mtext>cos</mtext><mi>β</mi></mrow></mtd><mtd><mrow><mo>−</mo><mtext>cos</mtext><mi>α</mi><mtext>sin</mtext><mi>β</mi></mrow></mtd></mtr><mtr><mtd><mrow><mtext>sin</mtext><mi>α</mi><mtext>cos</mtext><mi>β</mi></mrow></mtd><mtd><mrow><mtext>sin</mtext><mi>β</mi></mrow></mtd><mtd><mrow><mtext>cos</mtext><mi>α</mi><mtext>cos</mtext><mi>β</mi></mrow></mtd></mtr></mtable></mrow><mo>]</mo></mrow><mo>.</mo></mrow></math> The transformation results are shown in Fig. 2c 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: <math display="block" id="M11"><mrow><msup><mrow><mrow><mo>(</mo><mrow><msup><mi>x</mi><mo>′</mo></msup><mo>−</mo><msub><mi>x</mi><mi>s</mi></msub><msup><mrow/><mtext>"</mtext></msup></mrow><mo>)</mo></mrow></mrow><mn>2</mn></msup><mo>+</mo><msup><mrow><mrow><mo>(</mo><mrow><msup><mi>y</mi><mo>′</mo></msup><mo>−</mo><msub><mi>y</mi><mi>s</mi></msub><msup><mrow/><mtext>"</mtext></msup></mrow><mo>)</mo></mrow></mrow><mn>2</mn></msup><mo>=</mo><mfrac><mrow><mrow><mo>(</mo><mrow><mn>1</mn><mo>−</mo><msup><mrow><mtext>cos</mtext></mrow><mn>2</mn></msup><mi>θ</mi></mrow><mo>)</mo></mrow><msup><mrow><mrow><mo>(</mo><mrow><msup><mi>z</mi><mo>′</mo></msup><mo>−</mo><msub><mi>z</mi><mi>s</mi></msub><msup><mrow/><mtext>"</mtext></msup></mrow><mo>)</mo></mrow></mrow><mn>2</mn></msup></mrow><mrow><msup><mrow><mtext>cos</mtext></mrow><mn>2</mn></msup><mi>θ</mi></mrow></mfrac><mo>=</mo><msub><mi>R</mi><msup><mi>z</mi><mo>′</mo></msup></msub><msup><mrow/><mn>2</mn></msup><mo>.</mo></mrow></math> 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: <math display="block" id="M12"><mrow><mrow><mo>[</mo><mrow><mtable><mtr><mtd><mrow><msub><msup><mi>x</mi><mo>′</mo></msup><mi>s</mi></msub></mrow></mtd></mtr><mtr><mtd><mrow><msub><msup><mi>y</mi><mo>′</mo></msup><mi>s</mi></msub></mrow></mtd></mtr><mtr><mtd><mrow><msub><msup><mi>z</mi><mo>′</mo></msup><mi>s</mi></msub></mrow></mtd></mtr></mtable></mrow><mo>]</mo></mrow><mo>=</mo><mi>M</mi><mrow><mo>[</mo><mrow><mtable><mtr><mtd><mrow><msub><mi>x</mi><mi>s</mi></msub></mrow></mtd></mtr><mtr><mtd><mrow><msub><mi>y</mi><mi>s</mi></msub></mrow></mtd></mtr><mtr><mtd><mrow><msub><mi>z</mi><mi>s</mi></msub></mrow></mtd></mtr></mtable></mrow><mo>]</mo></mrow><mo>=</mo><mrow><mo>[</mo><mrow><mtable><mtr><mtd><mrow><msub><mi>x</mi><mi>s</mi></msub><mtext> cos </mtext><mi>α</mi><mo>−</mo><msub><mi>z</mi><mi>s</mi></msub><mtext> sin </mtext><mi>α</mi></mrow></mtd></mtr><mtr><mtd><mrow><mo>−</mo><msub><mi>x</mi><mi>s</mi></msub><mtext> sin </mtext><mi>α</mi><mtext> sin </mtext><mi>β</mi><mo>+</mo><msub><mi>y</mi><mi>s</mi></msub><mtext> cos </mtext><mi>β</mi><mo>−</mo><msub><mi>z</mi><mi>s</mi></msub><mtext>cos </mtext><mi>α</mi><mtext> sin </mtext><mi>β</mi></mrow></mtd></mtr><mtr><mtd><mrow><msub><mi>x</mi><mi>s</mi></msub><mtext> sin </mtext><mi>α</mi><mtext> cos </mtext><mi>β</mi><mo>+</mo><msub><mi>y</mi><mi>s</mi></msub><mtext>sin </mtext><mi>β</mi><mo>+</mo><msub><mi>z</mi><mi>s</mi></msub><mtext>cos </mtext><mi>α</mi><mtext> cos </mtext><mi>β</mi></mrow></mtd></mtr></mtable></mrow><mo>]</mo></mrow><mo>.</mo></mrow></math> 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 Fig. 2c . Meanwhile, <math id="M13"><mrow><mrow><mo>(</mo><mrow><msup><mi>x</mi><mo>′</mo></msup><mo>,</mo><mtext> </mtext><msup><mi>y</mi><mo>′</mo></msup><mo>,</mo><mtext> </mtext><msup><mi>z</mi><mo>′</mo></msup></mrow><mo>)</mo></mrow></mrow></math> are the coordinates on the conical surface. In the conical surface sampling process, <math id="M14"><msup><mi>z</mi><mo>′</mo></msup></math> can be determined in the direction of the axis of the standardized cone, and its parameter can be determined by uniform sampling as follows: <math display="block" id="M15"><mrow><msubsup><mi>Z</mi><mrow><mtext>plane</mtext></mrow><mn>2</mn></msubsup><mo>∼</mo><mi>U</mi><mrow><mo>(</mo><mrow><msubsup><mi>Z</mi><mrow><mtext>min</mtext></mrow><mn>2</mn></msubsup><mo>,</mo><mtext> </mtext><msubsup><mi>Z</mi><mrow><mtext>max</mtext></mrow><mn>2</mn></msubsup></mrow><mo>)</mo></mrow><mo>,</mo></mrow></math> where <math id="M16"><mrow><msub><mi>Z</mi><mrow><mtext>plane</mtext></mrow></msub></mrow></math> is the distance from the projection point of <math id="M17"><mrow><mrow><mo>(</mo><mrow><msup><mi>x</mi><mo>′</mo></msup><mo>,</mo><mtext> </mtext><msup><mi>y</mi><mo>′</mo></msup><mo>,</mo><mtext> </mtext><msup><mi>z</mi><mo>′</mo></msup></mrow><mo>)</mo></mrow></mrow></math> on the axis of the cone and the cone vertex, that is, the absolute value of <math id="M18"><mrow><msup><mi>z</mi><mo>′</mo></msup><mo>−</mo><msub><mi>z</mi><mi>s</mi></msub><msup><mrow/><mo>'</mo></msup></mrow></math> in Eq. (3); <math id="M19"><mrow><msub><mi>Z</mi><mrow><mtext>min</mtext></mrow></msub></mrow></math> and <math id="M20"><mrow><msub><mi>Z</mi><mrow><mtext>max</mtext></mrow></msub></mrow></math> are the minimum and maximum distances from the intersections of the sampling space and the axis of the cone to the cone vertex, respectively. Subsequently, according to Eq. (3), <math id="M21"><msup><mi>x</mi><mo>′</mo></msup></math> and <math id="M22"><msup><mi>y</mi><mo>′</mo></msup></math> can be determined using the following formulas: <math display="block" id="M23"><mrow><mtable><mtr><mtd><mrow><msup><mi>x</mi><mo>′</mo></msup><mo>=</mo><msub><msup><mi>x</mi><mo>′</mo></msup><mi>s</mi></msub><mo>+</mo><msub><mi>R</mi><msup><mi>z</mi><mo>′</mo></msup></msub><mi>cos</mi><mrow><mo>(</mo><mi>φ</mi><mo>)</mo></mrow></mrow></mtd></mtr><mtr><mtd><mrow><msup><mi>y</mi><mo>′</mo></msup><mo>=</mo><msub><msup><mi>y</mi><mo>′</mo></msup><mi>s</mi></msub><mo>+</mo><msub><mi>R</mi><msup><mi>z</mi><mo>′</mo></msup></msub><mi>sin</mi><mrow><mo>(</mo><mi>φ</mi><mo>)</mo></mrow></mrow></mtd></mtr></mtable><mtext>   </mtext><mrow><mo>(</mo><mrow><mtext> </mtext><mi>φ</mi><mo>∼</mo><mi>U</mi><mrow><mo>(</mo><mrow><mn>0</mn><mo>,</mo><mtext> </mtext><mn>2</mn><mi>π</mi></mrow><mo>)</mo></mrow></mrow><mo>)</mo></mrow><mo>.</mo></mrow></math> 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 <math id="M24"><mrow><mrow><mo>(</mo><mrow><msup><mi>x</mi><mo>′</mo></msup><mo>,</mo><mtext> </mtext><msup><mi>y</mi><mo>′</mo></msup><mo>,</mo><mtext> </mtext><msup><mi>z</mi><mo>′</mo></msup></mrow><mo>)</mo></mrow></mrow></math> is obtained, an inverse transformation can be applied to obtain the coordinates of <math id="M25"><mrow><mrow><mo>(</mo><mrow><mi>x</mi><mo>,</mo><mtext> y</mtext><mo>,</mo><mtext> z</mtext></mrow><mo>)</mo></mrow></mrow></math> in the laboratory coordinate system. Correction of the Compton scattering angle 3 The 3D-CZT CC prototype used in this study is fabricated by the Kromek group [ 44 44 , 45 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 <math id="M26"><mrow><msub><mi>δ</mi><mtext>E</mtext></msub></mrow></math> (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 Fig. 3a . The Compton scattering angle <math id="M28"><mi>θ</mi></math> is calculated using Eq. (7): <math display="block" id="M29"><mrow><mi>θ</mi><mo>=</mo><mtext>acos</mtext><mrow><mo>(</mo><mrow><mn>1</mn><mo>−</mo><mfrac><mrow><msub><mi>m</mi><mtext>e</mtext></msub><msup><mi>c</mi><mn>2</mn></msup><msub><mi>E</mi><mn>1</mn></msub></mrow><mrow><msub><mi>E</mi><mn>0</mn></msub><msub><mi>E</mi><mn>2</mn></msub></mrow></mfrac></mrow><mo>)</mo></mrow><mo>,</mo></mrow></math> where <math id="M30"><mrow><msub><mi>m</mi><mtext>e</mtext></msub><msup><mi>c</mi><mn>2</mn></msup></mrow></math> is the electron rest mass, <math id="M31"><mrow><msub><mi>E</mi><mn>0</mn></msub></mrow></math> is the incident gamma-ray energy, <math id="M32"><mrow><msub><mi>E</mi><mn>1</mn></msub></mrow></math> is the measured deposition energy of the recoil electrons by Compton scattering, and <math id="M33"><mrow><msub><mi>E</mi><mn>2</mn></msub></mrow></math> is 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 <math id="M34"><mi>θ</mi></math> [ 46 46 , 47 47 ]. The Gaussian blurring method for the cone angle can be found in Eq. (8) as follows: <math display="block" id="M35"><mrow><msup><mi>θ</mi><mo>′</mo></msup><mo>∼</mo><mi>N</mi><mrow><mo>(</mo><mrow><mi>θ</mi><mo>,</mo><mtext> </mtext><msub><mi>σ</mi><mi>θ</mi></msub></mrow><mo>)</mo></mrow><mo>,</mo></mrow></math> where <math id="M36"><mrow><msub><mi>σ</mi><mi>θ</mi></msub></mrow></math> represents the sigma of angular Gaussian blurring. For the conical surface sampling process described in the previous section, a new <math id="M37"><msup><mi>θ</mi><mo>′</mo></msup></math> is first generated according to Eq. (8) and subsequently introduced into Eqs. (3)–(6); this process is continuously repeated. Fig. 3b Fig. 3b presents the geometric relationship of Eq. (8), with <math id="M38"><mrow><msub><mi>σ</mi><mi>θ</mi></msub></mrow></math> individually calculated for each Compton event. This process guarantees a Gaussian distribution of the sampling points around the cone angle <math id="M39"><mi>θ</mi></math> . 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. Implementation of the MLEM iteration 2 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: <math display="block" id="M40"><mrow><msubsup><mi>f</mi><mi>j</mi><mrow><mi>l</mi><mo>+</mo><mn>1</mn></mrow></msubsup><mo>=</mo><mtext> </mtext><mfrac><mrow><msubsup><mi>f</mi><mi>j</mi><mi>l</mi></msubsup></mrow><mrow><msub><mi>S</mi><mi>j</mi></msub></mrow></mfrac><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>I</mi></munderover><mrow><mfrac><mrow><msub><mi>t</mi><mrow><mi>i</mi><mi>j</mi></mrow></msub></mrow><mrow><mstyle displaystyle="true"><msubsup><mo>∑</mo><mrow><mi>k</mi><mo>=</mo><mn>1</mn></mrow><mi>K</mi></msubsup><mrow><msubsup><mi>f</mi><mi>k</mi><mi>l</mi></msubsup><msub><mi>t</mi><mrow><mi>i</mi><mi>k</mi></mrow></msub></mrow></mstyle></mrow></mfrac></mrow></mstyle><mo>,</mo></mrow></math> where, <math id="M41"><mrow><mo> </mo><mi>j</mi></mrow></math> is the index of voxels in the imaging space, <math id="M42"><mrow><mo> </mo><mi>i</mi></mrow></math> is the index of detected Compton events, <math id="M43"><mrow><mo> </mo><msubsup><mi>f</mi><mi>j</mi><mi>l</mi></msubsup></mrow></math> is the intensity of the <math id="M44"><mrow><mi>j</mi><mtext>th</mtext></mrow></math> voxel in the imaging space after <math id="M45"><mi>l</mi></math> iterations, <math id="M46"><mrow><mo> </mo><msub><mi>S</mi><mi>j</mi></msub></mrow></math> is the relative probability of detecting the Compton events originating from the <math id="M47"><mrow><mi>j</mi><mtext>th</mtext></mrow></math> voxel, and the system matrix <math id="M48"><mrow><mo> </mo><msub><mi>t</mi><mrow><mi>i</mi><mi>j</mi></mrow></msub></mrow></math> corresponds to the probability of a photon emitted by the <math id="M49"><mrow><mi>j</mi><mtext>th</mtext></mrow></math> voxel to be detected as the <math id="M50"><mrow><mi>i</mi><mtext>th</mtext></mrow></math> Compton event. When l is equal to zero, the distribution of <math id="M51"><mrow><msubsup><mi>f</mi><mi>j</mi><mi>l</mi></msubsup></mrow></math> is the result of the Compton conical back-projection, that is, the precomputation of the LM-MLEM iteration. Moreover, <math id="M52"><mrow><msub><mi>S</mi><mi>j</mi></msub></mrow></math> was set to be uniform throughout the imaging space in this study [ 48 48 ]. The analytical calculation of the system matrix has been discussed by Maxim et al.; <math id="M53"><mrow><msub><mi>t</mi><mrow><mi>i</mi><mi>j</mi></mrow></msub></mrow></math> is finally positively related to the differential cross-section of the Compton scattering for the <math id="M54"><mrow><mi>i</mi><mi>t</mi><mi>h</mi></mrow></math> Compton event [ 34 34 ]. The forward projection of the LM-MLEM iteration contributes to the likelihood estimate of each Compton event: <math display="block" id="M55"><mrow><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>k</mi><mo>=</mo><mn>1</mn></mrow><mi>K</mi></munderover><mrow><msubsup><mi>f</mi><mi>k</mi><mi>l</mi></msubsup><msub><mi>t</mi><mrow><mi>i</mi><mi>k</mi></mrow></msub></mrow></mstyle><mo>.</mo></mrow></math> The forward projection of the LM-MLEM iteration is the summation of the system matrix element of each Compton event i multiplied by the image distribution. The <math id="M56"><mrow><msub><mi>t</mi><mrow><mi>i</mi><mi>k</mi></mrow></msub></mrow></math> depends only on index i and is nonzero on the set of voxels that contribute to the i th Compton event; thus, the following equations can be obtained: <math display="block" id="M57"><mrow><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>k</mi><mo>=</mo><mn>1</mn></mrow><mi>K</mi></munderover><mrow><msubsup><mi>f</mi><mi>k</mi><mi>l</mi></msubsup><msub><mi>t</mi><mrow><mi>i</mi><mi>k</mi></mrow></msub></mrow></mstyle><mo>=</mo><msub><mi>t</mi><mi>i</mi></msub><msubsup><mi>F</mi><mi>i</mi><mi>l</mi></msubsup><mo>,</mo></mrow></math> <math display="block" id="M58"><mrow><msubsup><mi>F</mi><mi>i</mi><mi>l</mi></msubsup><mo>=</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>k</mi><mo>=</mo><mn>1</mn></mrow><mi>K</mi></munderover><mrow><msubsup><mi>f</mi><mi>k</mi><mi>l</mi></msubsup></mrow></mstyle><mtext> </mtext><mrow><mo>(</mo><mrow><msub><mi>t</mi><mrow><mi>i</mi><mi>k</mi></mrow></msub><mo>≠</mo><mn>0</mn></mrow><mo>)</mo></mrow><mo>,</mo></mrow></math> where <math id="M59"><mrow><msub><mi>t</mi><mi>i</mi></msub></mrow></math> represents the factor positively related to the differential cross-section of the i th Compton event, and <math id="M60"><mrow><msubsup><mi>F</mi><mi>i</mi><mi>l</mi></msubsup></mrow></math> represents the sum of the intensities of the set of voxels that contribute to the i th Compton event. Based on the forward projection, the back-projection process of the LM-MLEM iteration can be expressed as follows: <math display="block" id="M61"><mrow><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>I</mi></munderover><mrow><mfrac><mrow><msub><mi>t</mi><mrow><mi>i</mi><mi>j</mi></mrow></msub></mrow><mrow><mstyle displaystyle="true"><msubsup><mo>∑</mo><mrow><mi>k</mi><mo>=</mo><mn>1</mn></mrow><mi>K</mi></msubsup><mrow><msubsup><mi>f</mi><mi>k</mi><mi>l</mi></msubsup><msub><mi>t</mi><mrow><mi>i</mi><mi>k</mi></mrow></msub></mrow></mstyle></mrow></mfrac></mrow></mstyle><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>I</mi></munderover><mo>=</mo></mstyle><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>I</mi></munderover><mrow><mfrac><mrow><msub><mi>t</mi><mrow><mi>i</mi><mi>j</mi></mrow></msub></mrow><mrow><msub><mi>t</mi><mi>i</mi></msub><msubsup><mi>F</mi><mi>i</mi><mi>l</mi></msubsup></mrow></mfrac></mrow></mstyle><mo>.</mo></mrow></math> The back-projection process of the LM-MLEM iteration is the summation of the ratio of the system matrix element of each Compton event i on voxel j to its forward projection, also called the back-projection factor of the <math id="M62"><mrow><mi>j</mi><mtext>th</mtext></mrow></math> voxel. In the calculation of the back-projection factor, <math id="M63"><mrow><msub><mi>t</mi><mrow><mi>i</mi><mi>j</mi></mrow></msub></mrow></math> represents the system matrix element for Compton event i , whose value is only related to index i ; therefore, the factor can be eliminated along with <math id="M64"><mrow><msub><mi>t</mi><mi>i</mi></msub></mrow></math> in the denominator. Subsequently, the back-projection factor can be calculated as follows: <math display="block" id="M65"><mrow><msubsup><mi>H</mi><mi>j</mi><mi>l</mi></msubsup><mo>=</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>I</mi></munderover><mrow><mfrac><mrow><msub><mi>t</mi><mrow><mi>i</mi><mi>j</mi></mrow></msub></mrow><mrow><msub><mi>t</mi><mi>i</mi></msub><msubsup><mi>F</mi><mi>i</mi><mi>l</mi></msubsup></mrow></mfrac></mrow></mstyle><mo>=</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>I</mi></munderover><mrow><mfrac><mn>1</mn><mrow><msubsup><mi>F</mi><mi>i</mi><mi>l</mi></msubsup></mrow></mfrac></mrow></mstyle><mtext> </mtext><mrow><mo>(</mo><mrow><msub><mi>t</mi><mi>i</mi></msub><mo>≠</mo><mn>0</mn></mrow><mo>)</mo></mrow><mo>.</mo></mrow></math> <math id="M66"><mrow><msubsup><mi>H</mi><mi>j</mi><mi>l</mi></msubsup></mrow></math> is the back-projection factor corresponding to the <math id="M67"><mrow><mi>j</mi><mtext>th</mtext></mrow></math> voxel in the <math id="M68"><mrow><mi>l</mi><mi>t</mi><mi>h</mi></mrow></math> LM-MLEM iteration. Considering that <math id="M69"><mrow><msub><mi>t</mi><mi>i</mi></msub></mrow></math> should be non-zero in the denominator, the computation of <math id="M70"><mrow><msubsup><mi>H</mi><mi>j</mi><mi>l</mi></msubsup></mrow></math> corresponds to the set of Compton events that contribute to the jth voxel. The LM-MLEM iterative form is finally expressed as follows: <math display="block" id="M71"><mrow><msubsup><mi>f</mi><mi>j</mi><mrow><mi>l</mi><mo>+</mo><mn>1</mn></mrow></msubsup><mo>=</mo><msubsup><mi>f</mi><mi>j</mi><mi>l</mi></msubsup><msubsup><mi>H</mi><mi>j</mi><mi>l</mi></msubsup><mo>.</mo></mrow></math> 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 × 10 5 . Parallel computing architecture for reconstruction 2 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 of <math id="M72"><mrow><msubsup><mi>f</mi><mi>j</mi><mn>0</mn></msubsup></mrow></math> . The total number of effective samplings of <math id="M73"><mrow><msub><mi>L</mi><mi>i</mi></msub></mrow></math> for the <math id="M74"><mrow><mi>i</mi><mtext>th</mtext></mrow></math> 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 length <math id="M75"><mi>L</mi></math> ( <math id="M76"><mrow><mi>L</mi><mo>=</mo><msub><mi>L</mi><mn>1</mn></msub><mo>+</mo><msub><mi>L</mi><mn>2</mn></msub><mo>+</mo><mo>⋯</mo><mo>+</mo><msub><mi>L</mi><mi>I</mi></msub></mrow></math> ) is constructed in the video memory to store the LM-MLEM iteration parameters, where <math id="M77"><mi>I</mi></math> is 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: For each event, <math id="M78"><mrow><msubsup><mi>F</mi><mi>i</mi><mi>l</mi></msubsup></mrow></math> in the forward projection of the LM-MLEM iteration is obtained by a thread on the GPU using the aforementioned parameter matrix. The back-projection factor <math id="M79"><mrow><msubsup><mi>H</mi><mi>j</mi><mi>l</mi></msubsup></mrow></math> is related to the set of Compton events corresponding to the jth voxel. Therefore, the fraction term <math id="M80"><mrow><mfrac><mn>1</mn><mrow><msubsup><mi>F</mi><mi>i</mi><mi>l</mi></msubsup></mrow></mfrac></mrow></math> in Equation (14) is accumulated in <math id="M81"><mrow><msubsup><mi>H</mi><mi>j</mi><mi>l</mi></msubsup></mrow></math> of each voxel in the voxel set corresponding to the i th event. Subsequently, when <math id="M82"><mrow><msubsup><mi>F</mi><mi>i</mi><mi>l</mi></msubsup></mrow></math> of all the events is calculated in parallel, <math id="M83"><mrow><msubsup><mi>H</mi><mi>j</mi><mi>l</mi></msubsup></mrow></math> of all the voxels are also accumulated, that is, the mapping of each voxel and its corresponding Compton event set is achieved. Finally, the intensity <math id="M84"><mrow><msubsup><mi>f</mi><mi>j</mi><mrow><mi>l</mi><mo>+</mo><mn>1</mn></mrow></msubsup></mrow></math> of each voxel in the imaging space after the iteration can be obtained by executing Equation (15). A simple description of the parallel architecture is presented in Fig. 4 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. Simulation and experimental settings 1 We performed both the Monte Carlo (MC) simulations and experiments to evaluate the performance of the proposed algorithm. Description of the Monte Carlo simulations 2 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 Fig. 2a , a multi-angle CC array composed of 3D-CZTs was constructed. Two 511 keV mono-energy point sources were simulated, as shown in Fig. 5a 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 Fig. 5b . The mouse phantom was referenced from a normal 28 g male mouse and generated by co-registered CT and cryosection images [ 49 49 ]. In the simulation, the phantom was divided into 380 × 992 × 208 voxels, and each voxel was 1 × 10 -3 mm 3 in 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 (C 12 H 18 O 4 , <math id="M85"><mi>ρ</mi></math> : 1.3 g/cm 3 ), 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. 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 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 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. Description of the mouse phantom experiment 2 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 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 multiangle 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 [ 18 F] sodium fluoride ([ 18 F]NaF), which was produced by JYAMS PET Research and Development Limited. The [ 18 F]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 50 ]. Results 1 Reconstruction results of the two-point source simulation 2 The reconstruction results for the two-point source were first evaluated using the FWHM along the three coordinate axes. Figure 7 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 the x -axis 22.16, 22.14, and 22.46 mm along the y -axis, and 22.42, 22.47, and 22.46 mm along the z -axis, respectively; the reconstruction times of these three methods are 16046.5, 1117.2, and 1223.3 s, respectively. Fig. 8 Fig. 8 presents the LM-MLEM 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 LM-MLEM is 7.22, 7.24, and 7.39 mm along the x -axis, 9.17, 9.24, and 9.42 mm along the y -axis, and 9.01, 8.93, and 9.19 mm along the z -axis, respectively; the reconstruction times of these three methods are 16339.7, 1357.9, and 1559.5 s, respectively. The reconstruction results shown in Figs. 7 Figs. 7 and 8 8 are the slices containing the radioactive source in the corresponding 3D reconstruction results. The calculation of the contrast ( V con ) is expressed by Eq. (16), <math display="block" id="M86"><mrow><msub><mi>V</mi><mrow><mtext>con</mtext></mrow></msub><mtext>=</mtext><mfrac><mrow><msub><mrow><mtext>Mean</mtext></mrow><mrow><mtext>ROI</mtext></mrow></msub><msub><mrow><mtext>-Mean</mtext></mrow><mrow><mtext>Background</mtext></mrow></msub></mrow><mrow><msub><mrow><mtext>Mean</mtext></mrow><mrow><mtext>Background</mtext></mrow></msub></mrow></mfrac><mo>,</mo></mrow></math> where <math id="M87"><mrow><msub><mrow><mtext>Mean</mtext></mrow><mrow><mtext>ROI</mtext></mrow></msub></mrow></math> is the average pixel intensity of the region of interest (ROI), and the <math id="M88"><mrow><msub><mrow><mtext>Mean</mtext></mrow><mrow><mtext>Background</mtext></mrow></msub></mrow></math> is 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 Fig. 8 , the calculated V con values 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 Table 1 . 1 1 Method 1 1 FWHM of x -axis (mm) 1 1 FWHM of y -axis (mm) 1 1 FWHM of z -axis (mm) 1 1 Reconstruction time (s) 1 1 SBP 1 1 15.53 1 1 22.16 1 1 22.42 1 1 16046.5 1 1 CSS-BP 1 1 15.19 1 1 22.14 1 1 22.47 1 1 1117.2 1 1 CSS-BP-SC 1 1 15.02 1 1 22.46 1 1 22.46 1 1 1223.3 1 1 SBP-based LM-MLEM(10 iterations) 1 1 7.22 1 1 9.17 1 1 9.01 1 1 16339.7 1 1 CSS-BP-based LM-MLEM(10 iterations) 1 1 7.24 1 1 9.24 1 1 8.93 1 1 1357.9 1 1 CSS-BP-SC-based LM-MLEM(10 iterations) 1 1 7.39 1 1 9.42 1 1 9.19 1 1 1537.0 Reconstruction results of the mouse phantom simulation 2 Figure 9 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 51 ]: <math display="block" id="M89"><mrow><msub><mrow><mtext>ARC</mtext></mrow><mi>i</mi></msub><mo>=</mo><mfrac><mrow><msub><mi>A</mi><mi>i</mi></msub><mstyle displaystyle="true"><msub><mo>∑</mo><mi>j</mi></msub><mrow><msub><mi>V</mi><mi>j</mi></msub></mrow></mstyle></mrow><mrow><msub><mi>V</mi><mi>i</mi></msub><msub><mi>A</mi><mi>T</mi></msub></mrow></mfrac><mo>,</mo></mrow></math> where <math id="M90"><mrow><msub><mi>A</mi><mi>i</mi></msub></mrow></math> and <math id="M91"><mrow><msub><mi>V</mi><mi>i</mi></msub></mrow></math> indicate the activity reconstructed in radioactive organ i and its volume, respectively, and <math id="M92"><mrow><msub><mi>A</mi><mtext>T</mtext></msub></mrow></math> is 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 GPU-accelerated 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 Table 2 lists the time consumption of the mouse phantom reconstruction process based on these three methods. 1 1 Method 1 1 Reconstruction time (s) 1 1 Relative reconstruction acceleration ratio 1 1 SBP-based LM-MLEM (10 iterations) 1 1 19771.2 1 1 1 1 1 CSS-BP-SC-based LM-MLEM (10 iterations) 1 1 1848.2 1 1 ~11 1 1 GPU-accelerated CSS-BP-SC-based LM-MLEM(10 iterations) 1 1 27.4 1 1 ~722 Reconstruction results of the mouse phantom experiment 2 Figure 10 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 time-consumptions 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 SBP-based, CSS-BP-SC-based, and GPU-accelerated CSS-BP-SC-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 Table 3 lists the time consumption of the mouse phantom reconstruction process based on these three methods. 1 1 Method 1 1 Reconstruction time (s) 1 1 Relative reconstruction acceleration ratio 1 1 SBP-based LM-MLEM (10 iterations) 1 1 10,750.7 1 1 1 1 1 CSS-BP-SC-based LM-MLEM (10 iterations) 1 1 1102.8 1 1 ~10 1 1 GPU-accelerated CSS-BP-SC-based LM-MLEM(10 iterations) 1 1 16.2 1 1 ~664 Discussion 1 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 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 Fig. 8 demonstrates that the CSS-BP-SC-based LM-MLEM improves the contrast V con by 9.3% compared to that of the CSS-BP-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. In this study, we designed a specialized parallel architecture for CSS-BP-SC-based LM-MLEM, as shown in Fig. 4 Fig. 4 . In the mouse phantom simulation, Fig. 9 Fig. 9 and Table 2 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 CSS-BP-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 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 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 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 8 , 9 9 , 10 10 , 11 11 , 12 12 , 23 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 53 , 54 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. Conclusion 1 In this study, we propose the use of a GPU-accelerated CSS-BP-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. References V. Filippou , C. Tsoumpas , Recent advances on the development of phantoms using 3D printing for imaging with CT, MRI, PET, SPECT, and ultrasound . Med. Phys. 45 , e740 - e760 ( 2018 ). doi: 10.1002/mp.13058 http://doi.org/10.1002/mp.13058 Recent advances on the development of phantoms using 3D printing for imaging with CT, MRI, PET, SPECT, and ultrasound C.R. Geng , Y. Ai , X.B. Tang et al. , Quantum dots enhanced Cerenkov luminescence imaging . Nucl. Sci. Tech. 30 , 71 ( 2019 ). doi: 10.1007/s41365-019-0599-x http://doi.org/10.1007/s41365-019-0599-x Quantum dots enhanced Cerenkov luminescence imaging X.W. Hu , D.D. Li , Y.J. Fu et al. , Advances in the application of radionuclide-labeled HER2 affibody for the diagnosis and treatment of ovarian cancer . Front. Oncol. 12 , 917439 ( 2022 ). doi: 10.3389/fonc.2022.917439 http://doi.org/10.3389/fonc.2022.917439 Advances in the application of radionuclide-labeled HER2 affibody for the diagnosis and treatment of ovarian cancer S. Motomura , Y. Kanayama , H. Haba et al. , Multiple molecular simultaneous imaging in a live mouse using semiconductor Compton camera . J. Anal. At. Spectrom 23 , 1089 ( 2008 ). doi: 10.1039/b802964d http://doi.org/10.1039/b802964d Multiple molecular simultaneous imaging in a live mouse using semiconductor Compton camera P. Kuchment and F. Terzioglu , Three-dimensional image reconstruction from compton camera data . SIAM. J. Imaging. Sci. 9 , 1708 - 1725 ( 2016 ). doi: 10.1137/16M107476X http://doi.org/10.1137/16M107476X Three-dimensional image reconstruction from compton camera data Y. Feng , A. Etxebeste , D. Sarrut et al. , 3-D reconstruction benchmark of a compton camera against a parallel-hole gamma camera on ideal data . IEEE Trans. Radiat. Plasma. Med. Sci. 4 , 479 - 488 ( 2020 ). doi: 10.1109/TRPMS.2019.2955745 http://doi.org/10.1109/TRPMS.2019.2955745 3-D reconstruction benchmark of a compton camera against a parallel-hole gamma camera on ideal data M. Fontana , D. Dauvergne , J. M. Létang et al. , Compton camera study for high efficiency SPECT and benchmark with Anger system . Phys. Med. Biol. 62 , 8794 - 8812 ( 2017 ). doi: 10.1088/1361-6560/aa926a http://doi.org/10.1088/1361-6560/aa926a Compton camera study for high efficiency SPECT and benchmark with Anger system Y. Suzuki , M. Yamaguchi , H. Odaka et al. , Three-dimensional and multienergy gamma-ray simultaneous imaging by using a Si/CdTe compton camera . Radiology 267 , 941 - 947 ( 2013 ). doi: 10.1148/radiol.13121194 http://doi.org/10.1148/radiol.13121194 Three-dimensional and multienergy gamma-ray simultaneous imaging by using a Si/CdTe compton camera M. Uenomach , M. Takahashi , K. Shimazoe et al. , Simultaneous in vivo imaging with PET and SPECT tracers using a Compton-PET hybrid camera . Sci. Rep. 11 , 1 - 11 ( 2021 ). doi: 10.1038/s41598-021-97302-7 http://doi.org/10.1038/s41598-021-97302-7 Simultaneous in vivo imaging with PET and SPECT tracers using a Compton-PET hybrid camera M. Sakai , M. Yamaguchi , Y. Nagao et al. , In vivo simultaneous imaging with 99m Tc and 18 F using a Compton camera . Phys. Med. Biol. 63 , 205006 ( 2018 ). doi: 10.1088/1361-6560/aae1d1 http://doi.org/10.1088/1361-6560/aae1d1 In vivo simultaneous imaging with 99mTc and 18F using a Compton camera K. Fujieda , J. Kataoka , S. Mochizuki et al. , First demonstration of portable Compton camera to visualize 223-Ra concentration for radionuclide therapy . Nucl. Instrum Methods Phys. Res. Sect A Accel. Spectrom. Detect. Assoc. Equip. 958 , 162802 ( 2020 ). doi: 10.1016/j.nima.2019.162802 http://doi.org/10.1016/j.nima.2019.162802 First demonstration of portable Compton camera to visualize 223-Ra concentration for radionuclide therapy T. Nakano , M. Sakai , K. Torikai et al. , Imaging of 99m Tc-DMSA and 18 F-FDG in humans using a Si/CdTe Compton camera . Phys. Med. Biol. 65 , 05LT01 ( 2020 ). doi: 10.1088/1361-6560/ab33d8 http://doi.org/10.1088/1361-6560/ab33d8 Imaging of 99m Tc-DMSA and 18 F-FDG in humans using a Si/CdTe Compton camera R. Todd , J. Nightingale , and D. Everett , A proposed γ camera . Nature 251 , 132 - 134 ( 1974 ). doi: 10.1038/251132a0 http://doi.org/10.1038/251132a0 A proposed γ camera M. Singh , R. R. Brechner , Experimental test-object study of electronically collimated SPECT . J. Nucl. Med. 31 , 178 - 186 ( 1990 ). Experimental test-object study of electronically collimated SPECT G. Llosa , G. Bernabeu , D. Burdette et al. , Development of a pre-clinical Compton probe prototype for prostate imaging . IEEE Symposium Conference Record Nuclear Science , 7 , 4168 - 4171 October 2004 Development of a pre-clinical Compton probe prototype for prostate imaging Z.Y. Yao , C.R. Shi , F. Tian et al. , Technical note: Rapid and high‐resolution deep learning–based radiopharmaceutical imaging with 3D‐CZT Compton camera and sparse projection data . Med. Phys. 49 , 7336 - 7346 ( 2022 ). doi: 10.1002/mp.15898 http://doi.org/10.1002/mp.15898 Technical note: Rapid and high‐resolution deep learning–based radiopharmaceutical imaging with 3D‐CZT Compton camera and sparse projection data F. Hueso-González , F. Fiedler , C. Golnik et al. , Compton camera and prompt gamma ray timing: Two methods for in vivo range assessment in proton therapy . Front. Oncol. 6 , 80 ( 2016 ). doi: 10.3389/fonc.2016.00080 http://doi.org/10.3389/fonc.2016.00080 Compton camera and prompt gamma ray timing: Two methods for in vivo range assessment in proton therapy D. Mackin , S. Peterson , S. Beddar et al. , Evaluation of a stochastic reconstruction algorithm for use in Compton camera imaging and beam range verification from secondary gamma emission during proton therapy . Phys. Med. Biol. 57 , 3537 - 3553 ( 2012 ). doi: 10.1088/0031-9155/57/11/3537 http://doi.org/10.1088/0031-9155/57/11/3537 Evaluation of a stochastic reconstruction algorithm for use in Compton camera imaging and beam range verification from secondary gamma emission during proton therapy R. K. Parajuli , M. Sakai , R. Parajuli et al. , Development and applications of compton camera—A review . Sensors 22 , 7374 ( 2022 ). doi: 10.3390/s22197374 http://doi.org/10.3390/s22197374 Development and applications of compton camera—A review Y. F. Yang , Y. Gono , S. Motomura et al. , A Compton camera for multitracer imaging . IEEE Trans. Nucl. Sci. 48 , 656 - 661 ( 2001 ). doi: 10.1109/23.940142 http://doi.org/10.1109/23.940142 A Compton camera for multitracer imaging B. Mehadji , M. Dupont , Y. Boursier et al. , Extension of the list-mode MLEM algorithm for poly-energetic imaging with a compton camera . IEEE Nuclear Science Symposium and Medical Imaging Conference Proceedings (NSS/MIC) , 1 - 8 November 2018 . doi: 10.1109/NSSMIC.2018.8824289 http://doi.org/10.1109/NSSMIC.2018.8824289 Extension of the list-mode MLEM algorithm for poly-energetic imaging with a compton camera Y.L. Liu , J.Q. Fu , Y.L. Li et al. , Preliminary results of a Compton camera based on a single 3D position-sensitive CZT detector . Nucl. Sci. Tech. 29 , 145 ( 2018 ). doi: 10.1007/s41365-018-0483-0 http://doi.org/10.1007/s41365-018-0483-0 Preliminary results of a Compton camera based on a single 3D position-sensitive CZT detector A. Kishimoto , J. Kataoka , T. Taya et al. , First demonstration of multi-color 3-D in vivo imaging using ultra-compact Compton camera . Sci. Rep. 7 , 2110 ( 2017 ). doi: 10.1038/s41598-017-02377-w http://doi.org/10.1038/s41598-017-02377-w First demonstration of multi-color 3-D in vivo imaging using ultra-compact Compton camera G. Rigaud and B. N. Hahn , Reconstruction algorithm for 3D Compton scattering imaging with incomplete data . Inverse. Probl. Sci. En. 29 , 967 - 989 ( 2021 ). doi: 10.1080/17415977.2020.1815723 http://doi.org/10.1080/17415977.2020.1815723 Reconstruction algorithm for 3D Compton scattering imaging with incomplete data Z.Y. Yao , Y.S Xiao , B. Wang et al. , Study of 3D fast Compton camera image reconstruction method by algebraic spatial sampling . Nucl. Instrum Methods Phys. Res. Sect A Accel. Spectrom. Detect. Assoc. Equip. 954 , 161345 ( 2020 ). doi: 10.1016/j.nima.2018.10.023 http://doi.org/10.1016/j.nima.2018.10.023 Study of 3D fast Compton camera image reconstruction method by algebraic spatial sampling S. J. Wilderman , W. L. Rogers , G. F. Knoll et al. , Fast algorithm for list mode back-projection of Compton scatter camera data . IEEE Trans. Nucl. Sci. 45 , 957 - 962 ( 1998 ). doi: 10.1109/23.682685 http://doi.org/10.1109/23.682685 Fast algorithm for list mode back-projection of Compton scatter camera data M. J. Cree and P. J. Bones , Towards direct reconstruction from a gamma camera based on Compton scattering . IEEE Trans. Med. Imaging. 13 , 398 - 407 ( 1994 ). doi: 10.1109/42.293932 http://doi.org/10.1109/42.293932 Towards direct reconstruction from a gamma camera based on Compton scattering R. Basko , G. L. Zeng , and G. T. Gullberg , Application of spherical harmonics to image reconstruction for the Compton camera . Phys. Med. Biol. 43 , 887 - 894 ( 1998 ). doi: 10.1088/0031-9155/43/4/016 http://doi.org/10.1088/0031-9155/43/4/016 Application of spherical harmonics to image reconstruction for the Compton camera M. Hirasawa and T. Tomitani , An analytical image reconstruction algorithm to compensate for scattering angle broadening in Compton cameras . Phys. Med. Biol. 48 , 1009 - 1026 ( 2003 ). doi: 10.1088/0031-9155/48/8/304 http://doi.org/10.1088/0031-9155/48/8/304 An analytical image reconstruction algorithm to compensate for scattering angle broadening in Compton cameras S. J. Wilderman , N. H. Clinthorne , J. A. Fessler et al. , List-mode maximum likelihood reconstruction of Compton scatter camera images in nuclear medicine . IEEE Nuclear Science Symposium and Medical Imaging Conference (Cat. No.98CH36255) , 3 , 1716 - 1720 November 1998 . doi: 10.1109/NSSMIC.1998.773871 http://doi.org/10.1109/NSSMIC.1998.773871 List-mode maximum likelihood reconstruction of Compton scatter camera images in nuclear medicine S. Takeda , H. Aono , S. Okuyama et al. , Experimental Results of the Gamma-Ray Imaging Capability With a Si/CdTe Semiconductor Compton Camera . IEEE Trans. Nucl. Sci. 56 , 783 - 790 ( 2009 ). doi: 10.1109/TNS.2008.2012059 http://doi.org/10.1109/TNS.2008.2012059 Experimental Results of the Gamma-Ray Imaging Capability With a Si/CdTe Semiconductor Compton Camera L. A. Shepp and Y. Vardi , Maximum Likelihood Reconstruction for Emission Tomography . IEEE Trans. Med. Imaging 1 , 113 - 122 ( 1982 ). doi: 10.1109/TMI.1982.4307558 http://doi.org/10.1109/TMI.1982.4307558 Maximum Likelihood Reconstruction for Emission Tomography H. H. Barrett , T. White , and L. C. Parra , List-mode likelihood . J. Opt. Soc. Am. A. 14 , 2914 ( 1997 ). doi: 10.1364/JOSAA.14.002914 http://doi.org/10.1364/JOSAA.14.002914 List-mode likelihood V. Maxim , X. Lojacono , E. Hilaire et al. , Probabilistic models and numerical calculation of system matrix and sensitivity in list-mode MLEM 3D reconstruction of Compton camera images . Phys. Med. Biol. 61 , 243 - 264 ( 2016 ). doi: 10.1088/0031-9155/61/1/243 http://doi.org/10.1088/0031-9155/61/1/243 Probabilistic models and numerical calculation of system matrix and sensitivity in list-mode MLEM 3D reconstruction of Compton camera images N. Kohlhase , T. Wegener , M. Schaar et al. , Capability of MLEM and OE to Detect Range Shifts With a Compton Camera in Particle Therapy . IEEE Trans. Radiat. Plasma Med. Sci. 4 , 233 - 242 ( 2020 ). doi: 10.1109/TRPMS.2019.2937675 http://doi.org/10.1109/TRPMS.2019.2937675 Capability of MLEM and OE to Detect Range Shifts With a Compton Camera in Particle Therapy Z.Y. Yao , Y.S. Xiao , Z.Q. Chen et al. , Compton-based prompt gamma imaging using ordered origin ensemble algorithm with resolution recovery in proton therapy . Sci. Rep. 9 , 1133 ( 2019 ). doi: 10.1038/s41598-018-37623-2 http://doi.org/10.1038/s41598-018-37623-2 Compton-based prompt gamma imaging using ordered origin ensemble algorithm with resolution recovery in proton therapy A. Andreyev , A. Celler , I. Ozsahin et al. , Resolution recovery for Compton camera using origin ensemble algorithm: Resolution recovery for Compton camera . Med. Phys. 43 , 4866 - 4876 ( 2016 ). doi: 10.1118/1.4959551 http://doi.org/10.1118/1.4959551 Resolution recovery for Compton camera using origin ensemble algorithm: Resolution recovery for Compton camera P. Després and X. Jia , A review of GPU-based medical image reconstruction . Physica Medica 42 , 76 - 92 ( 2017 ). doi: 10.1016/j.ejmp.2017.07.024 http://doi.org/10.1016/j.ejmp.2017.07.024 A review of GPU-based medical image reconstruction J.Y. Jiang , K. Li , S. Komarov et al. , Feasibility study of a point‐of‐care positron emission tomography system with interactive imaging capability . Med. Phys. 46 , 1798 - 1813 ( 2019 ). doi: 10.1002/mp.13397 http://doi.org/10.1002/mp.13397 Feasibility study of a point‐of‐care positron emission tomography system with interactive imaging capability X.Y Gu , L. Li , L. Wei et al. , Real-time reconstruction solution for positron emission mammography imaging-guided intervention . IEEE Nuclear Science Symposium and Medical Imaging Conference (NSS/MIC) , 1 - 5 November 2015 . doi: 10.1109/NSSMIC.2015.7582065 http://doi.org/10.1109/NSSMIC.2015.7582065 Real-time reconstruction solution for positron emission mammography imaging-guided intervention A.R. Zheng , Z.Y. Yao , and Y.S. Xiao , GPU Accelerated Stochastic Origin Ensemble Method With List-Mode Data for Compton Camera Imaging in Proton Therapy . IEEE Trans. Radiat. Plasma Med. Sci. 4 , 243 - 252 ( 2020 ). doi: 10.1109/TRPMS.2019.2929423 http://doi.org/10.1109/TRPMS.2019.2929423 GPU Accelerated Stochastic Origin Ensemble Method With List-Mode Data for Compton Camera Imaging in Proton Therapy VG. Nguyen and SJ. Lee , GPU-accelerated iterative reconstruction from Compton scattered data using a matched pair of conic projector and backprojector . Comput. Meth. Prog. Bio. 131 , 27 - 36 ( 2016 ). doi: 10.1016/j.cmpb.2016.04.012 http://doi.org/10.1016/j.cmpb.2016.04.012 GPU-accelerated iterative reconstruction from Compton scattered data using a matched pair of conic projector and backprojector Z.Y. Yao , Y.G. Yuan , J. Wu et al. , Rapid compton camera imaging for source terms investigation in the nuclear decommissioning with a subset-driven origin ensemble algorithm . Radiat. Phys. Chem. 197 , 110133 ( 2022 ). doi: 10.1016/j.radphyschem.2022.110133 http://doi.org/10.1016/j.radphyschem.2022.110133 Rapid compton camera imaging for source terms investigation in the nuclear decommissioning with a subset-driven origin ensemble algorithm F. Tian , C.R. Geng , Z.Y. Yao et al. , Radiopharmaceutical imaging based on 3D-CZT Compton camera with 3D-printed mouse phantom . Physica Medica 96 , 140 - 148 ( 2022 ). doi: 10.1016/j.ejmp.2022.03.005 http://doi.org/10.1016/j.ejmp.2022.03.005 Radiopharmaceutical imaging based on 3D-CZT Compton camera with 3D-printed mouse phantom Y. Li , P. Gong , X.B. Tang et al. , DOI correction for gamma ray energy reconstruction based on energy segment in 3D position-sensitive CdZnTe detectors . J. Inst. 17 , T03004 ( 2022 ). doi: 10.1088/1748-0221/17/03/T03004 http://doi.org/10.1088/1748-0221/17/03/T03004 DOI correction for gamma ray energy reconstruction based on energy segment in 3D position-sensitive CdZnTe detectors E. Yoshida , H. Tashima , K. Nagatsu et al. , Whole gamma imaging: a new concept of PET combined with Compton imaging . Phys. Med. Biol. 65 , 125013 ( 2020 ). doi: 10.1088/1361-6560/ab8e89 http://doi.org/10.1088/1361-6560/ab8e89 Whole gamma imaging: a new concept of PET combined with Compton imaging T. Ida , S. Motomura , M. Ueda et al. , Accurate modeling of event-by-event backprojection for a germanium semiconductor Compton camera for system response evaluation in the LM-ML-EM image reconstruction method . Jpn. J. Appl. Phys. 58 , 016002 ( 2019 ). doi: 10.7567/1347-4065/aae8e9 http://doi.org/10.7567/1347-4065/aae8e9 Accurate modeling of event-by-event backprojection for a germanium semiconductor Compton camera for system response evaluation in the LM-ML-EM image reconstruction method P. Solevi , E. Muñoz , C. Solaz et al. , Performance of MACACO Compton telescope for ion-beam therapy monitoring: first test with proton beams . Phys. Med. Biol. 61 , 5149 - 5165 ( 2016 ). doi: 10.1088/0031-9155/61/14/5149 http://doi.org/10.1088/0031-9155/61/14/5149 Performance of MACACO Compton telescope for ion-beam therapy monitoring: first test with proton beams B. Dogdas , D. Stout , A. F. Digimouse : a 3D whole body mouse atlas from CT and cryosection data . Phys. Med. Biol. 52 , 577 - 587 ( 2007 ). doi: 10.1088/0031-9155/52/3/003 http://doi.org/10.1088/0031-9155/52/3/003 : a 3D whole body mouse atlas from CT and cryosection data D. Shy and Z. He , Gamma-ray tracking for high energy gamma-ray imaging in pixelated CdZnTe . Nucl. Instrum Methods Phys. Res. Sect A Accel. Spectrom. Detect. Assoc. Equip. 954 , 161443 ( 2020 ). doi: 10.1016/j.nima.2018.10.121 http://doi.org/10.1016/j.nima.2018.10.121 Gamma-ray tracking for high energy gamma-ray imaging in pixelated CdZnTe E. Muñoz , A. Etxebeste , D. Dauvergne et al. , Imaging of polychromatic sources through Compton spectral reconstruction . Phys. Med. Biol. 67 , 195017 ( 2022 ). doi: 10.1088/1361-6560/ac92b9 http://doi.org/10.1088/1361-6560/ac92b9 Imaging of polychromatic sources through Compton spectral reconstruction Y. Kim , T. Lee , and W. Lee , Radiation measurement and imaging using 3D position sensitive pixelated CZT detector . Nucl. Eng. Technol. 51 , 1417 - 1427 ( 2019 ). doi: 10.1016/j.net.2019.03.009 http://doi.org/10.1016/j.net.2019.03.009 Radiation measurement and imaging using 3D position sensitive pixelated CZT detector H. Rohling , M. Priegnitz , S. Schoene et al. , Requirements for a Compton camera for in vivo range verification of proton therapy . Phys. Med. Biol. 62 , 2795 - 2811 ( 2017 ). doi: 10.1088/1361-6560/aa6068 http://doi.org/10.1088/1361-6560/aa6068 Requirements for a Compton camera for in vivo range verification of proton therapy C.H. Gong , X.B. Tang , D.Y. Shu et al. , Optimization of the Compton camera for measuring prompt gamma rays in boron neutron capture therapy . Appl. Radiat. Isot. 124 , 62 - 67 ( 2017 ). doi: 10.1016/j.apradiso.2017.03.014 http://doi.org/10.1016/j.apradiso.2017.03.014 Optimization of the Compton camera for measuring prompt gamma rays in boron neutron capture therapy All 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. 10.1007/s41365-023-01199-y This work was supported by the National Natural Science Foundation of China (No. 12220101005), Natural Science Foundation of Jiangsu Province (No. BK20220132), Primary Research and Development Plan of Jiangsu Province (No. BE2019002-3), Fundamental Research Funds for Central Universities (No. NG2022004), and the Foundation of the Graduate Innovation Center in NUAA (No. xcxjh20210613). 2022-11-20 2023-01-18 2023-02-08 2023-04-18 2023-06-07 © China Science Publishing & Media Ltd. (Science Press), Shanghai Institute of Applied Physics, the Chinese Academy of Sciences, ChineseNuclear Society and Springer Nature Singapore Pte Ltd. 2023 This is a post-peer-review, pre-copyedit version of an article published in Nuclear Science and Techniques. The final authenticated version is available online at: http://dx.doi.org/10.1007/s41365-023-01199-y http://dx.doi.org/10.1007/s41365-023-01199-y . Springer Nature or its licensor (e.g. a society or other partner) holds exclusive rights to this article under a publishing agreement with the author(s) or other rightsholder(s); author self-archiving of the accepted manuscript version of this article is solely governed by the terms of such publishing agreement and applicable law. 2023 38407957 The data that support the findings of this study are openly available in Science Data Bank at https://doi.org/10.57760/sciencedb.07763 and https://cstr.cn/31253.11.sciencedb.07763. Data availability Nuclear Science and Techniques 4 34