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 (3D-CZT 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.

Fig. 2Full-screen View

(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.

Equation (1) is the quadratic equation describing the Compton back-projection conical surface [23], where $\theta$ is the Compton scattering angle, $\overrightarrow{u}$ is the unit vector in the direction of the Compton axis, and $\left(x-x_s,y-y_s,z-z_s\right)$ represent the vectors connecting the points (x, y, z) on the conical surface and the Compton scattering points (x_s, y_s, z_s), $\begin{array}{l}{\left(\overrightarrow{u}\left(x-x_s,y-y_s,z-z_s\right)\right)}^2-(\left(x-x_s,y-y_s,z-z_s\right)\\ \cdot \left(x-x_s,y-y_s,z-z_s\right)){\text{cos}}^2\theta =0.\end{array}$ (1)

As shown in Fig. 2b, $\alpha$ and $\beta$ are the angles formed by the projection of $-\overrightarrow{u}$ in the XZ plane and the Z-axis, and the $-\overrightarrow{u}$ itself. The laboratory coordinate system can be transformed using the rotation matrix $M$ as follows: $M=\left[\begin{array}{ccc}\text{cos}\alpha & 0& -\text{sin}\alpha \\ -\text{sin}\alpha \text{sin}\beta & \text{cos}\beta & -\text{cos}\alpha \text{sin}\beta \\ \text{sin}\alpha \text{cos}\beta & \text{sin}\beta & \text{cos}\alpha \text{cos}\beta \end{array}\right].$ (2)

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: ${\left(x^{\prime }-x_s{}^{\text{"}}\right)}^2+{\left(y^{\prime }-y_s{}^{\text{"}}\right)}^2=\frac{\left(1-{\text{cos}}^2\theta \right){\left(z^{\prime }-z_s{}^{\text{"}}\right)}^2}{{\text{cos}}^2\theta }=R_{z^{\prime }}{}^2.$ (3)

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: $\left[\begin{array}{c}{x^{\prime }}_s\\ {y^{\prime }}_s\\ {z^{\prime }}_s\end{array}\right]=M\left[\begin{array}{c}{x}_s\\ y_s\\ z_s\end{array}\right]=\left[\begin{array}{c}{x}_s\text{ cos }\alpha -z_s\text{ sin }\alpha \\ -x_s\text{ sin }\alpha \text{ sin }\beta +y_s\text{ cos }\beta -z_s\text{cos }\alpha \text{ sin }\beta \\ x_s\text{ sin }\alpha \text{ cos }\beta +y_s\text{sin }\beta +z_s\text{cos }\alpha \text{ cos }\beta \end{array}\right].$ (4)

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, $\left(x^{\prime },y^{\prime },z^{\prime }\right)$ are the coordinates on the conical surface. In the conical surface sampling process, $z^{\prime }$ can be determined in the direction of the axis of the standardized cone, and its parameter can be determined by uniform sampling as follows: $Z_{\text{plane}}^2\sim U\left(Z_{\text{min}}^2,Z_{\text{max}}^2\right),$ (5) where $Z_{\text{plane}}$ is the distance from the projection point of $\left(x^{\prime },y^{\prime },z^{\prime }\right)$ on the axis of the cone and the cone vertex, that is, the absolute value of $z^{\prime }-z_s{}^{\text{'}}$ in Eq. (3); $Z_{\text{min}}$ and $Z_{\text{max}}$ 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), $x^{\prime }$ and $y^{\prime }$ can be determined using the following formulas: $\begin{array}{c}{x}^{\prime }={x^{\prime }}_s+R_{z^{\prime }}\cos \left(\phi \right)\\ y^{\prime }={y^{\prime }}_s+R_{z^{\prime }}\sin \left(\phi \right)\end{array}\left(\phi \sim U\left(0,2\pi \right)\right).$ (6)

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 $\left(x^{\prime },y^{\prime },z^{\prime }\right)$ is obtained, an inverse transformation can be applied to obtain the coordinates of $\left(x,\text{ y},\text{ z}\right)$ 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 ${\delta }_{\text{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.

Fig. 3Full-screen View

Correction of conical angle. (a) Energy resolution of the 3D-CZT CC. (b) Correction of the Compton scattering angle $\theta$.

The Compton scattering angle $\theta$ is calculated using Eq. (7): $\theta =\text{acos}\left(1-\frac{m_{\text{e}}{c}^2{E}_1}{E_0{E}_2}\right),$ (7) where $m_{\text{e}}{c}^2$ is the electron rest mass, $E_0$ is the incident gamma-ray energy, $E_1$ is the measured deposition energy of the recoil electrons by Compton scattering, and $E_2$ 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 $\theta$ [46, 47]. The Gaussian blurring method for the cone angle can be found in Eq. (8) as follows: ${\theta }^{\prime }\sim N\left(\theta ,{\sigma }_{\theta }\right),$ (8) where ${\sigma }_{\theta }$ represents the sigma of angular Gaussian blurring. For the conical surface sampling process described in the previous section, a new ${\theta }^{\prime }$ is first generated according to Eq. (8) and subsequently introduced into Eqs. (3)–(6); this process is continuously repeated. Fig. 3b presents the geometric relationship of Eq. (8), with ${\sigma }_{\theta }$ individually calculated for each Compton event. This process guarantees a Gaussian distribution of the sampling points around the cone angle $\theta$.

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

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: $f_j^{l+1}=\frac{f_j^l}{S_j}{\displaystyle \sum _{i=1}^I\frac{t_{ij}}{{\displaystyle {\sum }_{k=1}^K{f}_k^l{t}_{ik}}}},$ (9) where, $j$ is the index of voxels in the imaging space, $i$ is the index of detected Compton events, $f_j^l$ is the intensity of the $j\text{th}$ voxel in the imaging space after $l$ iterations, $S_j$ is the relative probability of detecting the Compton events originating from the $j\text{th}$ voxel, and the system matrix $t_{ij}$ corresponds to the probability of a photon emitted by the $j\text{th}$ voxel to be detected as the $i\text{th}$ Compton event.

When l is equal to zero, the distribution of $f_j^l$ is the result of the Compton conical back-projection, that is, the precomputation of the LM-MLEM iteration. Moreover, $S_j$ was 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.; $t_{ij}$ is finally positively related to the differential cross-section of the Compton scattering for the $ith$ Compton event [34].

The forward projection of the LM-MLEM iteration contributes to the likelihood estimate of each Compton event: ${\displaystyle \sum _{k=1}^K{f}_k^l{t}_{ik}}.$ (10)

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 $t_{ik}$ depends only on index i and is nonzero on the set of voxels that contribute to the ith Compton event; thus, the following equations can be obtained: ${\displaystyle \sum _{k=1}^K{f}_k^l{t}_{ik}}=t_i{F}_i^l,$ (11) $F_i^l={\displaystyle \sum _{k=1}^K{f}_k^l}\left(t_{ik}\ne 0\right),$ (12) where $t_i$ represents the factor positively related to the differential cross-section of the ith Compton event, and $F_i^l$ represents the sum of the intensities of the set of voxels that contribute to the ith Compton event. Based on the forward projection, the back-projection process of the LM-MLEM iteration can be expressed as follows: ${\displaystyle \sum _{i=1}^I\frac{t_{ij}}{{\displaystyle {\sum }_{k=1}^K{f}_k^l{t}_{ik}}}}{\displaystyle \sum _{i=1}^I=}{\displaystyle \sum _{i=1}^I\frac{t_{ij}}{t_i{F}_i^l}}.$ (13)

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 $j\text{th}$ voxel. In the calculation of the back-projection factor, $t_{ij}$ 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 $t_i$ in the denominator. Subsequently, the back-projection factor can be calculated as follows: $H_j^l={\displaystyle \sum _{i=1}^I\frac{t_{ij}}{t_i{F}_i^l}}={\displaystyle \sum _{i=1}^I\frac{1}{F_i^l}}\left(t_i\ne 0\right).$ (14)

$H_j^l$ is the back-projection factor corresponding to the $j\text{th}$ voxel in the $lth$ LM-MLEM iteration. Considering that $t_i$ should be non-zero in the denominator, the computation of $H_j^l$ corresponds to the set of Compton events that contribute to the jth voxel. The LM-MLEM iterative form is finally expressed as follows: $f_j^{l+1}=f_j^l{H}_j^l.$ (15)

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⁵.

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 of $f_j^0$. The total number of effective samplings of $L_i$ for the $i\text{th}$ 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 $L$ ($L=L_1+L_2+\cdots +L_I$) is constructed in the video memory to store the LM-MLEM iteration parameters, where $I$ 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, $F_i^l$ 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 $H_j^l$ is related to the set of Compton events corresponding to the jth voxel. Therefore, the fraction term $\frac{1}{F_i^l}$ in Equation (14) is accumulated in $H_j^l$ of each voxel in the voxel set corresponding to the ith event. Subsequently, when $F_i^l$ of all the events is calculated in parallel, $H_j^l$ 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 $f_j^{l+1}$ 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. 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.

Fig. 4Full-screen View

Architecture of parallel computing

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 multi-angle 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^-3 mm³ 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₁₂H₁₈O₄, $\rho$: 1.3 g/cm³), 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. 5Full-screen View

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

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.

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 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 [¹⁸F] sodium fluoride ([¹⁸F]NaF), which was produced by JYAMS PET Research and Development Limited. The [¹⁸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].

Fig. 6Full-screen View

(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.

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 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 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 and 8 are the slices containing the radioactive source in the corresponding 3D reconstruction results.

Fig. 7Full-screen View

(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.

Fig. 8Full-screen View

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

The calculation of the contrast (V_con) is expressed by Eq. (16), $V_{\text{con}}\text{=}\frac{{\text{Mean}}_{\text{ROI}}{\text{-Mean}}_{\text{Background}}}{{\text{Mean}}_{\text{Background}}},$ (16) where ${\text{Mean}}_{\text{ROI}}$ is the average pixel intensity of the region of interest (ROI), and the ${\text{Mean}}_{\text{Background}}$ 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, 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.

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]: ${\text{ARC}}_i=\frac{A_i{\displaystyle {\sum }_j{V}_j}}{V_i{A}_T},$ (17) where $A_i$ and $V_i$ indicate the activity reconstructed in radioactive organ i and its volume, respectively, and $A_{\text{T}}$ 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 lists the time consumption of the mouse phantom reconstruction process based on these three methods.

Fig. 9Full-screen View

(Color online) Reconstruction results of the SBP-based, CSS-BP-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.

Reconstruction results of the mouse phantom experiment

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 lists the time consumption of the mouse phantom reconstruction process based on these three methods.

Fig. 10Full-screen View

(Color online) Reconstruction results of the SBP-based, CSS-BP-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.