A method and system for improving reconstruction speed and signal-to-noise ratio of a cascaded gamma imaging system

Abstract

The application provides a method and system for improving the reconstruction speed and signal-to-noise ratio of a cascade gamma imaging system, and relates to the technical field of nuclear medical imaging. The method comprises the following steps: according to the acquired slit coincidence event data, combining the gamma photon interaction information and the collimator type to encode the detector module, and dividing the coincidence event data into multiple subsets according to the energy information and the encoding; directly back-projecting the coincidence events in each subset to determine the projection intersection points, and mapping the projection intersection points into corresponding voxel numbers; defining a region of interest with the voxel position represented by the voxel number as the center, and determining the boundary of the region of interest; according to the region of interest voxels, calculating the system transmission matrix corresponding to each pair of slit coincidence events in each data subset, based on the region of interest voxels, the divided data subsets and the system transmission matrix, using the ordered subset expectation maximization algorithm to perform image iterative reconstruction; thereby improving the reconstruction speed and signal-to-noise ratio of the cascade gamma imaging system.

Description

Technical Field

[0001] This invention relates to the field of nuclear medicine imaging technology, and in particular to a method and system for improving the reconstruction speed and signal-to-noise ratio of a cascaded gamma imaging system. Background Technology

[0002] Therapeutic integration is a medical model that combines diagnosis and treatment. It utilizes a single targeted molecular carrier to simultaneously load diagnostic and therapeutic radionuclides, thereby achieving synergistic effects of precise tumor imaging and targeted therapy. For example... 177 Lu、 131 I, 111 In、 47 Sc、 161 Nuclides such as Tb can release both beta and gamma rays during their decay. Gamma rays are suitable for single-photon emission computed tomography (SPECT), and can be used for diagnostic tracing and drug distribution studies; while beta rays carry a single negative charge and are easily absorbed by surrounding tissues, thus having potential applications in tumor therapy. Among many nuclides, 177 Lu, due to its moderate half-life (approximately 6.7 days) and ease of production and transportation, has been widely used in clinical and preclinical research. This is thanks to its dual role as a gamma emitter (primarily with energies of 208 keV and 113 keV) and a medium-energy beta emitter (maximum energy 0.497 MeV, penetration distance in tissue approximately 2 mm), combining SPECT imaging capabilities with effective tumor cell killing while causing minimal damage to surrounding normal tissues. This achieves a "dual-purpose" approach, making it considered an ideal therapeutic radionuclide. Currently, 177 Lu primarily binds to prostate-specific membrane antigen (PSMA) ligands or somatostatin analogues (DOTATATE) and is used in the diagnosis and treatment of prostate cancer and neuroendocrine tumors.

[0003] However, traditional SPECT systems suffer from low sensitivity and limited spatial resolution, making it difficult to achieve real-time dynamic imaging. Furthermore, 177 During the decay of Lu, there is approximately a 3.2% probability of generating cascaded gamma photon pairs that are highly correlated in time and space. Compared to traditional single-photon imaging, this phenomenon provides more imaging information and helps improve the image signal-to-noise ratio. To address this, researchers have proposed a spatiotemporal coincidence imaging method based on slit collimation, aiming to improve the image signal-to-noise ratio. 177Imaging detection is performed using cascaded gamma photons generated by Lu decay. During this imaging process, coincident gamma photon pairs passing through the aperture collimator and the slit collimator are simultaneously recorded, i.e., aperture-slit coincidence events. Based on these coincidence events, back projections can be performed along the center point of the narrowest part of the aperture collimator and the center line of the narrowest part of the slit collimator, respectively, to establish back projection lines and back projection planes. Then, by calculating the coordinates of their intersection points, the aperture collimator can be roughly determined. 177 The decay location of Lu. This type of direct backprojection reconstruction method has high real-time performance, but poor signal-to-noise ratio and relies on a large amount of aperture matching data.

[0004] In addition, some studies have employed the Maximum Likelihood Expectation Maximization (MLEM) algorithm to iteratively reconstruct data acquired by cascaded systems. This method can achieve a good signal-to-noise ratio with relatively small amounts of data. However, similar to positron emission tomography (PET) and SPECT, there are two ways to calculate the system transfer matrix in cascaded coincidence imaging: one is to pre-calculate and store the entire system matrix offline and call it directly during iteration. Although this method has a fast reconstruction speed, it requires loading the huge system matrix into the runtime memory, which consumes significant memory resources and is therefore difficult to implement in practical applications. The other method is to calculate the system response matrix for each coincidence event in real time during each iteration. Although this method avoids the storage of the system matrix, it is slow and inefficient. Furthermore, in PET and SPECT, the system matrix calculation only needs to consider the intersection of the response line and the voxels in the imaging field of view; while in cascaded coincidence imaging, each coincidence event needs to traverse all voxels in the field of view and calculate the corresponding response value, which further reduces its calculation speed, efficiency, and the signal-to-noise ratio of the reconstructed image.

[0005] Therefore, how to improve the computational efficiency and speed of the system matrix in the cascaded gamma-photon coincidence imaging system (hereinafter referred to as the cascaded gamma imaging system) and accelerate the convergence process of iterative reconstruction has become a key problem that urgently needs to be solved. Summary of the Invention

[0006] To address the technical problem of slow reconstruction speed and time-consuming system transfer matrix calculations due to the need to traverse the entire field of view voxels in existing cascaded gamma imaging iterative reconstruction algorithms, this invention provides a method and system for improving the reconstruction speed and signal-to-noise ratio of cascaded gamma imaging systems. The technical solution is as follows:

[0007] On the one hand, a method for improving the reconstruction speed and signal-to-noise ratio of a cascaded gamma imaging system is provided, including: S1. Based on the collected coincidence event data, the detector module is encoded in combination with the gamma photon interaction information and the collimator type, and the coincidence event data is divided into multiple subsets according to the energy information and the encoding; among them, the coincidence event is the aperture coincidence event; S2. Perform direct back projection on the coincidence events in each subset to determine the projection intersection points and map them to the corresponding voxel numbers; S3. Delineate the region of interest centered on the voxel position represented by the voxel number, and determine the boundary of the region of interest; S4. Calculate the system transfer matrix corresponding to each pair of aperture coincidence events in each data subset. Based on the partitioned data subsets and the system transfer matrix, perform iterative image reconstruction using the ordered subset expectation maximization algorithm. When calculating the system transfer matrix, only the response probability between the coincidence event and each voxel in its corresponding region of interest is calculated. During the iterative reconstruction process, when calculating the projection and back projection values ​​of each pair of coincidence events, only the voxels in the region of interest are traversed.

[0008] Furthermore, the step of encoding the detector module based on the collected coincidence event data, combined with gamma photon interaction information and collimator type, and dividing the coincidence event data into multiple subsets based on energy information and encoding includes: S11. Based on the position of the interaction between the γ photon and the crystal and the type of collimator through which the γ photon passes, the position of each detector module or detector unit is encoded. S12. Based on the energy information of the γ photons in the aperture coincidence event and the position code, divide the aperture coincidence event data into N subsets according to a preset rule.

[0009] Furthermore, the step of directly back-projecting the coincident events in each subset to determine the projection intersection points and mapping them to the corresponding voxel numbers includes: S21. Traverse all hole-slit coincidence event pairs in each subset, perform direct backprojection to establish a backprojection model, and calculate the coordinates of the intersection point of the corresponding backprojection line and the backprojection plane in the field of view space. S22. Based on the coordinates of the intersection point of the aperture matching event pair in the field of view space, map them to the corresponding field of view voxel number.

[0010] Furthermore, the process of traversing all aperture-slit coincidence event pairs in each subset, performing direct backprojection to establish a backprojection model, and calculating the coordinates of the intersection point of the corresponding backprojection line and the backprojection plane in the field of view includes: S211. Based on the interaction position of the γ photon passing through the aperture collimator with the crystal in the aperture coincidence event, calculate the coordinates of the center point of the corresponding crystal unit, and establish the straight line equation of the back projection line pointing from the center point of the crystal unit to the center point of the aperture collimator. S212. Determine the vector of the centerline of the slit collimator; S213. Based on the vector of the center line of the slit collimator, the coordinates of the crystal unit center, and the coordinates of the center point of the slit collimator, establish the plane equation of the back projection surface; S214. Based on the plane equation of the back-projection surface and the straight line equation of the back-projection line, solve for the three-dimensional coordinates of their intersection point. .

[0011] Further, mapping the intersection coordinates of the aperture coincidence event pairs in the field of view space to the corresponding field of view voxel numbers includes: S221. Pre-define the field of view to be reconstructed and divide it into a regular grid in three-dimensional space. Each grid cell is called a voxel. Let the FOV be... X The boundary of the direction is The number of voxels is ; Y Directional boundary is The number of voxels is The boundary in the Z direction is The number of voxels is The dimensions of each voxel in the three dimensions are: ; In the formula, , , They are respectively the voxels along the field of view X , Y , Z Dimensions in three directions; S222, Calculate the intersection point The voxel number corresponding to the intersection point is obtained by dividing the offset relative to the lower boundary of the FOV by the voxel size and rounding down. ; In the formula, This represents the floor function. , , They are respectively X , Y , Z The voxel index in the direction, and satisfying , , ; S223. For intersections that fall outside the FOV boundary, the corresponding events are removed.

[0012] Furthermore, the step of defining the region of interest centered on the voxel position represented by the voxel number and determining the boundary of the region of interest includes: S31. Taking the voxel position represented by the voxel number as the center, define a preset size within its neighborhood. m × n × l The voxel region is taken as the region of interest, where, m , n , l respectively along X axis, Y axis, Z Number of voxels along the axial direction; S32. Determine the boundary of the region of interest in the overall field of view voxel space based on the voxel space of the overall imaging field of view, the voxel number corresponding to the intersection point, and the size of the region of interest.

[0013] Further, determining the boundary of the region of interest in the overall field of view voxel space based on the voxel space of the overall imaging field of view, the voxel number corresponding to the intersection point, and the size of the region of interest includes: Assume the entire imaging field of view voxel space is in X axis, Y axis, Z The voxel index ranges along the axis are as follows: X axis: ,common Individual factors; Y axis: ,common Individual factors; Z axis: ,common Individual factors; Calculate the number of voxels extending outwards from the central voxel index in each dimension, and define the half-width in each direction as: ; In the formula, , , The regions of interest are respectively in X , Y , Z Half width in the direction; The theoretical voxel boundary of the region of interest is then: ; ; In the formula, , , The regions of interest are respectively in X , Y , Z Theoretical starting index in the direction; , , The regions of interest are respectively in X , Y , Z Theoretical termination index in direction, for Central voxel numbering; The obtained theoretical voxel boundaries are corrected; the correction formula is as follows: ; ; After correction, the region of interest is in X Axis effective voxel index range: Region of interest is Y Axis effective voxel index range: Region of interest is Z Axis effective voxel index range: .

[0014] Furthermore, the ordered subset expectation-maximization algorithm is used for iterative image reconstruction, with the following formula: ; In the formula, Voxel values ​​for the field of view and Global numbering of field voxels. For the number of iterations, It is the event number. The total number of events, The total number of voxels in the field of view. For the first The cascaded gamma photons produced by nuclear decay in the voxel of the field of view are the first The probability of an event. For the first The cascaded gamma photons produced by nuclear decay in the voxel of the field of view are the first The probability of an event. For the first A subset of data, superscript For the first In the first iteration Data subsets Iteration, It is the projected value. , , Indicates the region of interest.

[0015] On the other hand, a system for improving the reconstruction speed and signal-to-noise ratio of a cascaded gamma imaging system is provided. This system is used to implement the aforementioned method for improving the reconstruction speed and signal-to-noise ratio of a cascaded gamma imaging system. The system includes: The subset partitioning module is used to encode the detector module based on the collected coincidence event data, combined with the gamma photon interaction information and collimator type, and to divide the coincidence event data into multiple subsets based on the energy information and the encoding. The voxel mapping module is used to perform direct back projection on the coincidence events in each subset to determine the projection intersection points and map them to the corresponding voxel numbers; The region determination module is used to delineate the region of interest centered on the voxel position represented by the voxel number, and to determine the boundary of the region of interest; The iterative reconstruction module is used to calculate the system transfer matrix corresponding to each pair of aperture coincidence events in each data subset. Based on the partitioned data subsets and the system transfer matrix, the ordered subset expectation maximization algorithm is used for iterative image reconstruction. When calculating the system transfer matrix, only the response probability between the coincidence event and each voxel in the corresponding region of interest is calculated. During the iterative reconstruction process, when calculating the projection and back projection values ​​of each pair of coincidence events, only the voxels in the region of interest are traversed.

[0016] The beneficial effects of the technical solutions provided by the embodiments of the present invention include at least the following: (1) By classifying the aperture coincidence events according to the collimator type they pass through, the encoding of the interacting detector modules, and the energy value of the γ photons, the aperture coincidence event data is divided into several subsets. The ordered subset expectation maximization algorithm is used to iteratively reconstruct each subset, which significantly accelerates the iteration convergence and thus effectively improves the reconstruction speed of the cascaded γ imaging system.

[0017] (2) During the iterative reconstruction process, only voxels in the region of interest are traversed when calculating the projection and backprojection values ​​of each pair of coincident events. This helps to limit the possible location range of radioactive source decay, effectively narrowing the computational range by limiting the region of interest. Thus, it solves the problems of traditional iterative reconstruction algorithms, such as the need to traverse all voxels in the field of view to calculate the system transfer matrix, high time consumption of projection and backprojection operations, and slow convergence speed. This further improves the computational speed of iterative reconstruction and enhances the signal-to-noise ratio of the system imaging and the overall image signal-to-noise ratio.

[0018] (3) By directly backprojecting the slit coincidence event to solve for the intersection point, and establishing the region of interest centered on the intersection point, only the response probability values ​​of the event and the voxels within the region of interest need to be calculated when calculating the system transfer matrix elements. In this way, it is avoided that all field-of-view voxels need to be traversed when calculating the system transfer matrix elements, thereby significantly improving the computational efficiency of the system transfer matrix and helping to further improve the reconstruction speed of the cascaded gamma imaging system. Attached Figure Description

[0019] To more clearly illustrate the technical solutions in the embodiments of the present invention, the accompanying drawings used in the description of the embodiments will be briefly introduced below. Obviously, the accompanying drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0020] Figure 1 This is a flowchart of a method for improving the reconstruction speed and signal-to-noise ratio of a cascaded gamma imaging system according to an embodiment of the present invention; Figure 2 This is a schematic diagram of the position numbering of the cascaded gamma imaging system provided in an embodiment of the present invention; Figure 3 This is a flowchart of subset partitioning in an embodiment of the method for improving the reconstruction speed and signal-to-noise ratio of a cascaded gamma imaging system provided by the present invention; Figure 4 This is a schematic diagram of dividing data subsets according to the energy information of γ photons in the cascaded γ imaging system provided in the embodiment of the present invention; Figure 5 This is a schematic diagram of dividing data subsets based on the aperture collimator in the cascaded gamma imaging system provided in this embodiment of the invention; Figure 6 This is a schematic diagram of dividing data subsets based on the slit collimator in the cascaded gamma imaging system provided in the embodiments of the present invention; Figure 7 This is a schematic diagram of a cascaded gamma imaging system provided in an embodiment of the present invention, which uses the aperture collimator as the standard and considers energy information to divide the data subset; Figure 8 This is a schematic diagram of a cascaded gamma imaging system provided in an embodiment of the present invention, which uses a slit collimator as the standard and considers energy information to divide the data subset; Figure 9 This is a schematic diagram of the back projection and voxel mapping process provided in an embodiment of the present invention; Figure 10 This is a schematic diagram of the region of interest determined in the cascaded gamma imaging system provided in the embodiments of the present invention; Figure 11 This is a schematic diagram of the process for determining the region of interest provided in an embodiment of the present invention; Figure 12 This is a schematic diagram of the iterative reconstruction process provided in an embodiment of the present invention; Figure 13 This is a system block diagram of the cascaded gamma imaging system provided in the embodiments of the present invention, which improves the reconstruction speed and signal-to-noise ratio. Figure reference numerals: 1-10 are the detector numbers on the ring; 11, the system's field of view; 12, the slit collimator; 13, the aperture collimator; 14, the gamma photon detected by the detector corresponding to the aperture collimator; 15, the gamma photon detected by the detector corresponding to the slit collimator; 16, the detector module number corresponding to the aperture collimator on the ring; 17, the detector module number corresponding to the slit collimator on the ring; 18, the detector module number on the ring corresponding to the aperture collimator for low-energy gamma photons; 19, the detector module number on the ring corresponding to the slit collimator for high-energy gamma photons; 20, the detector module number on the ring corresponding to the aperture collimator for high-energy gamma photons; 21, the detector module number on the ring corresponding to the slit collimator for low-energy gamma photons; 22, the back projection line; 23, the back projection plane; 24, the intersection of the back projection line and the back projection plane; 25, the region of interest S. Detailed Implementation

[0021] The technical solution of the present invention will now be described with reference to the accompanying drawings.

[0022] In embodiments of the present invention, words such as "exemplarily," "for example," etc., are used to indicate that something is an example, illustration, or description. Any embodiment or design described as "exemplary" in the present invention should not be construed as being more preferred or advantageous than other embodiments or designs. Specifically, the use of the word "exemplary" is intended to present the concept in a concrete manner. Furthermore, in embodiments of the present invention, the meaning expressed by "and / or" can be both, or either one.

[0023] In the embodiments of this invention, the terms "image" and "picture" may sometimes be used interchangeably. It should be noted that, without emphasizing the distinction between them, they convey the same meaning. Similarly, the terms "of," "corresponding (relevant)," and "corresponding" may sometimes be used interchangeably. It should be noted that, without emphasizing the distinction between them, they convey the same meaning.

[0024] In this embodiment of the invention, sometimes a subscript such as W1 may be written in a non-subscript form such as W1. When the difference is not emphasized, the meaning they express is the same.

[0025] To make the technical problems, technical solutions and advantages of the present invention clearer, a detailed description will be given below in conjunction with the accompanying drawings and specific embodiments.

[0026] like Figure 1As shown, this embodiment of the invention provides a method for improving the reconstruction speed and signal-to-noise ratio of a cascaded gamma imaging system. Essentially, it is a reconstruction method based on subset partitioning and the expected value maximization of ordered subsets of the region of interest. The processing flow of this method may include the following steps: S1. Based on the collected coincidence event data, the detector module is encoded in combination with the gamma photon interaction information and the collimator type, and the coincidence event data is divided into multiple subsets according to the energy information and the encoding. To better understand this embodiment, a brief description of the cascaded gamma imaging system will be provided first, such as... Figure 2 As shown, the structure of the cascaded gamma imaging system mainly includes: detector modules numbered 1-10, a field of view 11, a slit collimator 12, an aperture collimator 13, and a data acquisition module, etc.; among which, The total number of detector modules is 40, arranged in a ring structure with 4 axes, with 10 detector modules evenly distributed on each ring; in practical applications, the position and number of detector modules can be arranged according to actual needs, and this embodiment does not limit the number and arrangement of detector modules. The detector module measures 54.9 mm × 54.9 mm × 10 mm, and the detector material is cerium-doped gadolinium aluminum gallium garnet (Ce, GAGG:Ce). In practical applications, the size of the detector module and the type of detector material can be set as needed. This embodiment does not limit the size of the detector module or the type of material. The field of view is a cube with dimensions of 40 mm × 40 mm × 40 mm. In practical applications, the shape and size of the field of view can be set as needed. This embodiment does not limit the selection of the shape and size of the field of view. In this embodiment, tungsten can be selected as the collimator material; in practical applications, the materials of the slit collimator and the hole collimator can be any one of high-density metals such as tungsten and lead. This embodiment does not limit the selection of materials for the slit collimator and the hole collimator. The slit collimators and hole collimators must correspond one-to-one with the detector modules on the ring, but the ratio of the number of hole collimators to slit collimators can be any value. This embodiment does not limit the ratio of the number of hole collimators to slit collimators. In this embodiment, the ratio of the number of hole collimators to slit collimators is 1:1.

[0027] The data acquisition module described above can use the PETsys system. However, this embodiment does not limit the specific type of data acquisition module.

[0028] like Figure 3 As shown, step S1 may specifically include the following steps: S11. Based on the position of the interaction between the γ photon and the crystal (referring to the detector module) and the type of collimator through which the γ photon passes, the position is encoded for each detector module or detector unit. In this embodiment, the location of the interaction between the γ photon and the crystal specifically includes the x, y, and z axis coordinates of the interaction between the γ photon and the detector, the detector module number of the ring where the interaction location is located, the axial detector module number, and the crystal array number in each detector module.

[0029] Furthermore, based on the coordinate information of the interaction between the γ photon and the detector, the module number on the ring, the axial module number, the crystal array number in each detector module, and the collimator type through which the γ photon passes, each detector module or detector unit is position-coded according to a preset rule.

[0030] In this embodiment, the position encoding of the detector module or detector unit can be performed based on the following methods: on-ring detector module, axial detector module, all detector modules, crystal array in each detector module, crystal array in on-ring detector module, crystal array in axial detector module, or crystal array in all detector modules.

[0031] In this embodiment, each detector module on the ring can be encoded (modules with the same axis have the same encoding). Of course, this embodiment does not limit the specific method of position encoding. The rules for encoding the position of each detector module on the ring include both clockwise and counterclockwise encoding. This embodiment uses clockwise encoding for each detector module on the ring. This embodiment does not limit the selection of encoding rules for each detector module on the ring.

[0032] In this embodiment, when encoding the clockwise positions of each detector module on the ring, any detector module on the ring can be set as the encoding starting point. In this embodiment, detector module 1 at the highest point on the ring is used as the encoding starting point. Of course, this embodiment does not limit the selection of the encoding starting point.

[0033] S12. Based on the energy information of the γ photons in the aperture coincidence event and the position code, divide the aperture coincidence event data into N subsets according to a preset rule.

[0034] Since the cascaded gamma imaging system involved in this embodiment uses a slit collimator to detect cascaded gamma photons, the approximate location of the nuclide decay for each pair of slit coincidence events can be obtained through direct back projection, thus ensuring sampling completeness for each pair of slit coincidence events. Therefore, the slit coincidence event data can be divided into multiple subsets, and an ordered subset expectation-maximization algorithm can be used to accelerate the reconstruction convergence speed. Based on this, according to the energy information and position encoding of the gamma photon pairs in the slit coincidence event data, the slit coincidence event data is divided into N subsets according to a pre-set rule.

[0035] In this embodiment, the preset rule may include: relying solely on the energy information of γ photons (such as...) Figure 4 As shown), using the aperture collimator as the standard (e.g.) Figure 5 As shown), using the seam straightener as the standard (e.g.) Figure 6 As shown), using the aperture collimator as the standard and considering energy information (such as... Figure 7 As shown), using the slit collimator as the standard and considering energy information (such as... Figure 8 The data subsets are divided using various methods, including using the aperture collimator as the standard and a one-to-one correspondence, using the slit collimator as the standard and a one-to-one correspondence, using the aperture collimator as the standard and a one-to-one correspondence while considering energy information, and using the slit collimator as the standard and a one-to-one correspondence while considering energy information. The number of data subsets obtained by different division methods may also differ. This embodiment uses the method of using the aperture collimator as the standard and considering energy information to divide the aperture-slit matching data for reconstruction into 10 subsets. It is understood that this embodiment does not limit the specific choice of data subset division rules.

[0036] In this embodiment, the principle of data subset partitioning based on the aperture collimator and considering energy information is as follows: According to the on-ring detector encoding, γ-photon events in aperture-slit coincidence events that pass through the same aperture collimator are divided into the same subset, resulting in 5 subsets; based on this, the energy information of the γ-photon is further considered, i.e., the energy of the γ-photon passing through the aperture collimator may be 113 keV or 208 keV, thus expanding the total number of data subsets to 10. Specifically, as follows... Figure 7 As shown, Figure 7 In the diagram, 18 represents the detector module number on the ring corresponding to the low-energy γ photon through-aperture collimator, 19 represents the detector module number on the ring corresponding to the high-energy γ photon through-slit collimator, 20 represents the detector module number on the ring corresponding to the high-energy γ photon through-aperture collimator, and 21 represents the detector module number on the ring corresponding to the low-energy γ photon through-slit collimator.

[0037] S2. Perform direct backprojection on the coincident events in each subset to determine the projection intersection points, and map them to the corresponding voxel numbers; such as Figure 9 As shown, the specific steps may include: S21. Traverse all hole-slit coincidence event pairs in each subset, perform direct backprojection to establish a backprojection model, and calculate the coordinates of the intersection point of the corresponding backprojection line and the backprojection plane in the field of view space; specifically, this may include the following steps: S211. Based on the interaction position of the γ photon through the aperture collimator 13 with the crystal in the aperture coincidence event, calculate the coordinates of the center point of the corresponding crystal unit, and establish the straight line equation of the back projection line 22 pointing from the center point of the crystal unit to the center point of the aperture collimator. S212, Determine the vector of the centerline of the slit collimator 12; S213. Based on the vector of the center line of the slit collimator, the coordinates of the crystal unit center, and the coordinates of the center point of the slit collimator, establish the plane equation of the back projection plane 23. S214. Based on the plane equation of the back-projection surface and the straight line equation of the back-projection line, solve for the three-dimensional coordinates of their intersection point 24. .

[0038] S22. Map the coordinates of the intersection points of the aperture matching event pairs in the field of view space to the corresponding field of view voxel numbers; specifically, this may include the following steps: S221. Pre-define the field of view to be reconstructed and divide it into a regular grid in three-dimensional space. Each grid cell is called a voxel. Let the field of view (FOV) be at... X The boundary of the direction is The number of voxels is ; Y Directional boundary is The number of voxels is The boundary in the Z direction is The number of voxels is The dimensions of each voxel in the three dimensions are: ; In the formula, , , They are respectively the voxels along the field of view X , Y , Z Dimensions in three directions; In this embodiment, the field of view in the cascaded gamma imaging system is along... X , Y , Z Number of voxels in direction , , It can be set according to actual needs. For example, the number of voxels. , , All are set to 50.

[0039] S222, Calculate the intersection point The voxel number corresponding to the intersection point is obtained by dividing the offset relative to the lower boundary of the FOV by the voxel size and rounding down. ; In the formula, This represents the floor function. , , They are respectively X , Y , Z The voxel index in the direction, and satisfying , , ; S223. For intersections that fall outside the FOV boundary, the corresponding events are removed.

[0040] S3, such as Figure 10 As shown, a region of interest 25 is defined centered on the voxel position represented by the voxel number, and the boundary of the region of interest is determined; as... Figure 11 As shown, the specific steps may include: S31. The voxel number corresponding to the intersection point obtained from step S22. Centered on the voxel position represented by the voxel number, a preset size is defined within its neighborhood. m × n × l The voxel region is taken as the region of interest, where, m , n , l respectively along X axis, Y axis, Z Number of voxels along the axial direction; In this embodiment, the preset size m × n × l The number of voxel regions can be set according to actual imaging needs, and can be odd or even. When m , n , l average When the number is odd, the central voxel of the region of interest is... And it is symmetrically distributed in all three dimensions; when m , n , l When there is an even number of voxels, the center point is located between adjacent voxels. In this embodiment... m , n , l All values ​​are set to 5. Of course, this embodiment does not limit the selection of the size of the region of interest.

[0041] S32. Based on the voxel space of the overall imaging field of view and the voxel number corresponding to the intersection point. and the size of the region of interest m × n × l Determine the boundary of the region of interest in the overall field of view voxel space, that is, its boundary in... X , Y , Z The starting and ending voxel indices are located in three dimensions; specifically, this may include the following steps: Assume the entire imaging field of view voxel space is in X axis, Y axis, Z The voxel index ranges along the axis are as follows: X axis: ,common Individual factors; Y axis: ,common Individual factors; Z axis: ,common Individual factors; Calculate the number of voxels extending outwards from the central voxel index in each dimension, and define the half-width in each direction as: ; In the formula, , , The regions of interest are respectively in X , Y , Z Half width in the direction; The theoretical voxel boundary of the region of interest is then: ; ; In the formula, , , The regions of interest are respectively in X , Y , Z Theoretical starting index in the direction; , , The regions of interest are respectively in X , Y , Z Theoretical termination index in direction, for Central voxel numbering; Since the calculated theoretical voxel boundaries may exceed the actual index range of the entire imaging field of view voxel space, boundary checks and corrections are necessary to ensure the validity of all voxel indices. The correction formula is as follows: ; ; After correction, the region of interest is in X Axis effective voxel index range: Region of interest is Y Axis effective voxel index range: Region of interest is Z Axis effective voxel index range: .

[0042] S4. Calculate the system transfer matrix corresponding to each pair of aperture coincidence events in each data subset. Based on the partitioned data subsets and the system transfer matrix, perform iterative image reconstruction using the ordered subset expectation-maximization algorithm. Specifically, when calculating the system transfer matrix, only the response probability between the coincidence event and each voxel within its corresponding region of interest is calculated. During the iterative reconstruction process, when calculating the projection and back-projection values ​​for each pair of coincidence events, only the voxels within the region of interest are traversed. Figure 12 As shown, the specific steps may include: S41. Calculate the system transfer matrix corresponding to each pair of aperture coincidence events in each data subset. When calculating the system transfer matrix elements corresponding to each pair of aperture coincidence events in each subset, only the elements corresponding to that coincidence event and of size [missing element] are traversed. m × n × l Voxels within the region of interest are used to calculate the response probability value associated with the coincidence event; In this embodiment, for each pair of aperture coincidence events in all subsets, the corresponding three-dimensional region of interest is determined, and the size of this region is [size missing]. m × n × l Voxels. When calculating the system transfer matrix, it is only necessary to traverse the voxels within the corresponding region of interest and calculate the response probability of the coincidence event with each voxel within the region of interest, without having to calculate its probability with all voxels in the entire field of view. This method reduces the computational complexity of elements in the system transfer matrix from O(n log n) to O(n log n). Reduced to O( m × n × l ),in This represents the total number of voxels in the entire image, significantly reducing computation time and improving computational efficiency.

[0043] In this embodiment, the response probability value can be calculated based on methods such as point spread function models, geometric factors, Monte Carlo simulations, and manual measurements. This embodiment uses a method based on geometric factors for calculation; however, other methods can also be used, and there is no limitation.

[0044] S42, based on the division N Using a subset of data and the calculated system transfer matrix, an ordered subset expectation-maximization algorithm is employed for iterative image reconstruction. In each iteration, when calculating the projection value for each pair of coincident events and performing the back-projection operation, only the data of the specified size are traversed according to the region of interest. m × n × l Voxels within the region of interest.

[0045] In this embodiment, based on the N data subsets and their corresponding system transmission matrices (probability values ​​obtained through real-time calculation), an ordered subset expectation-maximization algorithm is used for iterative image reconstruction. The formula is as follows: ; In the formula, Voxel values ​​for the field of view and Global numbering of field voxels. For the number of iterations, It is the event number. The total number of events, The total number of voxels in the field of view. For the first The cascaded gamma photons produced by nuclear decay in the voxel of the field of view are the first The probability of an event. For the first The cascaded gamma photons produced by nuclear decay in the voxel of the field of view are the first The probability of an event. For the first A subset of data, superscript For the first In the first iteration Data subsets Iteration, It is the projected value. and They are the first In the first iteration Data subsets After the first iteration and before the second iteration j voxel values ​​of the field of view Indicates the first In the first iteration Data subsets voxel values ​​of the field of view before iteration , , Indicates the region of interest.

[0046] In this embodiment, the image needs to be initialized before using the ordered subset expectation maximization algorithm for iterative image reconstruction.

[0047] Specifically, the initialization and reconstruction of the image can be achieved by setting it to a completely uniform image, based on the results of direct backprojection, or by randomization. This embodiment uses the method based on the results of direct backprojection to initialize the image; however, this embodiment does not limit the specific method of initializing and reconstructing the image.

[0048] In each iteration, the iterative algorithm reconstructs each subset sequentially according to a preset subset order, specifically including the following steps: For the subset currently being processed, a forward projection is first performed based on the current image estimation. Calculate the forward projection estimate ; Subsequently, by using the actual projected values With the corresponding projection estimate Compare and the ratio Multiplying by the probability yields a reverse projection value; where, , This is the ratio of the actual projected value to the corresponding estimated projected value. Next, all backprojection values ​​in this subset are summed and backprojected to the image space to update the image with respect to the entire field of view voxels; After processing one subset, move on to the next subset, and continue until all subsets have been traversed, thus completing one full iteration. Finally, determine whether the preset number of iterations has been reached or the convergence criterion has been met. If not, continue to the next round of iterations.

[0049] In this embodiment, during the iteration of the aforementioned data subset, when calculating the forward and back projection operations for each coincident event, local operations are performed based on the region of interest. Specifically:

[0050] Calculate the forward projection for a given coincident event. At that time, only the region of interest corresponding to the event is accumulated. m × n × l Current image estimates of each voxel within the image Its probability value The product of, i.e. , Indicates the region of interest; In the back projection operation, the projection ratio is also... The product of the projection ratio and the system probability matrix elements is limited to the region of interest. Specifically, it involves multiplying the projection ratio by the system probability value of the corresponding voxel within the region. .

[0051] In this embodiment, the calculation of projection and backprojection is limited to the region of interest during the iteration process, thereby avoiding calculations on the entire image and reducing the computational complexity from O(n log n) at the global scale. Reduced to the local scale O( m × n × l This significantly improves the overall efficiency of image reconstruction.

[0052] This embodiment, while maintaining the quality of the reconstructed image, limits the calculation scope of the system matrix calculation, projection, and backprojection operations from the global image to the region of interest, significantly reducing memory requirements and computational complexity, and achieving efficient and fast image reconstruction.

[0053] In this embodiment, the convergence criterion for the above iteration can be set based on indicators such as contrast-to-noise ratio (CNR), structural similarity index measure (SSIM), normalized root mean square error (NRMSE), peak signal-to-noise ratio (PSNR), contrast recovery coefficient (CRC), and mean squared error (MSE). This embodiment uses both the structural similarity index and the contrast recovery coefficient to determine whether the iteration converges. However, this embodiment does not limit the method for determining the iteration convergence criterion.

[0054] In this embodiment, a parallel computing strategy can also be adopted during the iteration process to perform parallel processing of system matrix calculation, forward projection and back projection operations for different coincidence events. This is especially suitable for graphics processing units (GPUs) or distributed computing environments, and can greatly accelerate the reconstruction process.

[0055] In summary, the beneficial effects of the method for improving the reconstruction speed and signal-to-noise ratio of the cascaded gamma imaging system provided in this embodiment of the invention include at least the following: (1) By classifying the aperture coincidence events according to the collimator type they pass through, the encoding of the interacting detector modules, and the energy value of the γ photons, the aperture coincidence event data is divided into several subsets. The ordered subset expectation maximization algorithm is used to iteratively reconstruct each subset, which significantly accelerates the iteration convergence and thus effectively improves the reconstruction speed of the cascaded γ imaging system.

[0056] (2) During the iterative reconstruction process, only voxels in the region of interest are traversed when calculating the projection and backprojection values ​​of each pair of coincident events. This helps to limit the possible location range of radioactive source decay, effectively narrowing the computational range by limiting the region of interest. Thus, it solves the problems of traditional iterative reconstruction algorithms, such as the need to traverse all voxels in the field of view to calculate the system transfer matrix, high time consumption of projection and backprojection operations, and slow convergence speed. This further improves the computational speed of iterative reconstruction and enhances the signal-to-noise ratio of the system imaging and the overall image signal-to-noise ratio.

[0057] (3) By directly backprojecting the slit coincidence event to solve for the intersection point, and establishing the region of interest centered on the intersection point, only the response probability values ​​of the event and the voxels within the region of interest need to be calculated when calculating the system transfer matrix elements. In this way, it is avoided that all field-of-view voxels need to be traversed when calculating the system transfer matrix elements, thereby significantly improving the computational efficiency of the system transfer matrix and helping to further improve the reconstruction speed of the cascaded gamma imaging system.

[0058] Figure 13 This is a system block diagram illustrating an exemplary embodiment for improving the reconstruction speed and signal-to-noise ratio of a cascaded gamma imaging system. The system is used to implement a method for improving the reconstruction speed and signal-to-noise ratio of a cascaded gamma imaging system. (Refer to...) Figure 13 The system includes a subset partitioning module 10, a voxel mapping module 20, a region determination module 30, and an iterative reconstruction module 40. Among them:

[0059] The subset partitioning module 10 is used to encode the detector module based on the collected coincidence event data, combined with the gamma photon interaction information and the collimator type, and to divide the coincidence event data into multiple subsets based on the energy information and the encoding. The voxel mapping module 20 is used to perform direct back projection on the coincidence events in each subset to determine the projection intersection point and map it to the corresponding voxel number. The region determination module 30 is used to delineate the region of interest centered on the voxel position represented by the voxel number, and to determine the boundary of the region of interest; The iterative reconstruction module 40 is used to calculate the system transfer matrix corresponding to each pair of aperture coincidence events in each data subset. Based on the partitioned data subsets and the system transfer matrix, the ordered subset expectation maximization algorithm is used for image iterative reconstruction. When calculating the system transfer matrix, only the response probability between the coincidence event and each voxel in the corresponding region of interest is calculated. During the iterative reconstruction process, when calculating the projection and back projection values ​​of each pair of coincidence events, only the voxels in the region of interest are traversed.

[0060] The above embodiments can be implemented, in whole or in part, by software, hardware (such as circuits), firmware, or any other combination thereof. When implemented using software, the above embodiments can be implemented, in whole or in part, as a computer program product. The computer program product includes one or more computer instructions or computer programs. When the computer instructions or computer programs are loaded or executed on a computer, all or part of the processes or functions described in the embodiments of the present invention are generated. The computer can be a general-purpose computer, a special-purpose computer, a computer network, or other programmable device. The computer instructions can be stored in a computer-readable storage medium or transmitted from one computer-readable storage medium to another. For example, the computer instructions can be transmitted from one website, computer, server, or data center to another website, computer, server, or data center via wired (e.g., infrared, wireless, microwave, etc.) means. The computer-readable storage medium can be any available medium that a computer can access or a data storage device such as a server or data center that includes one or more sets of available media. The available medium can be a magnetic medium (e.g., floppy disk, hard disk, magnetic tape), an optical medium (e.g., DVD), or a semiconductor medium. A semiconductor medium can be a solid-state drive.

[0061] It should be understood that the term "and / or" in this article is merely a description of the relationship between related objects, indicating that three relationships can exist. For example, A and / or B can represent: A alone, A and B simultaneously, or B alone. A and B can be singular or plural. Additionally, the character " / " in this article generally indicates an "or" relationship between the preceding and following related objects, but it can also represent an "and / or" relationship. Please refer to the context for a more accurate understanding.

[0062] In this invention, "at least one" means one or more, and "more than one" means two or more. "At least one of the following" or similar expressions refer to any combination of these items, including any combination of a single item or a plurality of items. For example, at least one of a, b, or c can represent: a, b, c, ab, ac, bc, or abc, where a, b, and c can be single or multiple.

[0063] It should be understood that, in various embodiments of the present invention, the order of the above-mentioned process numbers does not imply the order of execution. The execution order of each process should be determined by its function and internal logic, and should not constitute any limitation on the implementation process of the embodiments of the present invention.

[0064] Those skilled in the art will recognize that the units and algorithm steps of the various examples described in conjunction with the embodiments disclosed herein can be implemented in electronic hardware, or a combination of computer software and electronic hardware. Whether these functions are implemented in hardware or software depends on the specific application and design constraints of the technical solution. Those skilled in the art can use different methods to implement the described functions for each specific application, but such implementations should not be considered beyond the scope of this invention.

[0065] Those skilled in the art will clearly understand that, for the sake of convenience and brevity, the specific working processes of the devices, apparatuses, and units described above can be referred to the corresponding processes in the foregoing method embodiments, and will not be repeated here.

[0066] In the several embodiments provided by this invention, it should be understood that the disclosed devices, apparatuses, and methods can be implemented in other ways. For example, the apparatus embodiments described above are merely illustrative; for instance, the division of units is only a logical functional division, and in actual implementation, there may be other division methods. For example, multiple units or components may be combined or integrated into another device, or some features may be ignored or not executed. Furthermore, the coupling or direct coupling or communication connection shown or discussed may be through some interfaces; the indirect coupling or communication connection between devices or units may be electrical, mechanical, or other forms.

[0067] The units described as separate components may or may not be physically separate. The components shown as units may or may not be physical units; that is, they may be located in one place or distributed across multiple network units. Some or all of the units can be selected to achieve the purpose of this embodiment according to actual needs.

[0068] In addition, the functional units in the various embodiments of the present invention can be integrated into one processing unit, or each unit can exist physically separately, or two or more units can be integrated into one unit.

[0069] If the aforementioned functions are implemented as software functional units and sold or used as independent products, they can be stored in a computer-readable storage medium. Based on this understanding, the technical solution of this invention, or the part that contributes to the prior art, or a part of the technical solution, can be embodied in the form of a software product. This computer software product is stored in a storage medium and includes several instructions to cause a computer device (which may be a personal computer, server, or network device, etc.) to execute all or part of the steps of the methods described in the various embodiments of this invention. The aforementioned storage medium includes various media capable of storing program code, such as USB flash drives, portable hard drives, read-only memory (ROM), random access memory (RAM), magnetic disks, or optical disks.

[0070] The above description is merely a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the technical scope disclosed in the present invention should be included within the scope of protection of the present invention. Therefore, the scope of protection of the present invention should be determined by the scope of the claims.

Claims

1. A method for improving the reconstruction speed and signal-to-noise ratio of a cascaded gamma imaging system, characterized in that, The method includes: S1. Based on the collected coincidence event data, the detector module is encoded in combination with the gamma photon interaction information and the collimator type, and the coincidence event data is divided into multiple subsets according to the energy information and the encoding; among them, the coincidence event is the aperture coincidence event; S2. Perform direct back projection on the coincidence events in each subset to determine the projection intersection points and map them to the corresponding voxel numbers; S3. Delineate the region of interest centered on the voxel position represented by the voxel number, and determine the boundary of the region of interest; S4. Calculate the system transfer matrix corresponding to each pair of aperture coincidence events in each data subset. Based on the partitioned data subsets and the system transfer matrix, perform iterative image reconstruction using the ordered subset expectation maximization algorithm. When calculating the system transfer matrix, only the response probability between the coincidence event and each voxel in its corresponding region of interest is calculated. During the iterative reconstruction process, when calculating the projection and back projection values ​​of each pair of coincidence events, only the voxels in the region of interest are traversed.

2. The method for improving the reconstruction speed and signal-to-noise ratio of a cascaded gamma imaging system according to claim 1, characterized in that, The process involves encoding the detector module based on the collected coincidence event data, combined with gamma photon interaction information and collimator type, and dividing the coincidence event data into multiple subsets based on energy information and encoding, including: S11. Based on the position of the interaction between the γ photon and the crystal and the type of collimator through which the γ photon passes, the position of each detector module or detector unit is encoded. S12. Based on the energy information of the γ photons in the aperture coincidence event and the position code, divide the aperture coincidence event data into N subsets according to a preset rule.

3. The method for improving the reconstruction speed and signal-to-noise ratio of a cascaded gamma imaging system according to claim 1, characterized in that, The step of directly back-projecting the coincident events in each subset to determine the projection intersection points and mapping them to the corresponding voxel numbers includes: S21. Traverse all hole-slit coincidence event pairs in each subset, perform direct backprojection to establish a backprojection model, and calculate the coordinates of the intersection point of the corresponding backprojection line and the backprojection plane in the field of view space. S22. Based on the coordinates of the intersection point of the aperture matching event pair in the field of view space, map them to the corresponding field of view voxel number.

4. The method for improving the reconstruction speed and signal-to-noise ratio of a cascaded gamma imaging system according to claim 3, characterized in that, The process of traversing all aperture-slit coincidence event pairs in each subset, performing direct backprojection to establish a backprojection model, and calculating the coordinates of the intersection point of the corresponding backprojection line and the backprojection plane in the field of view includes: S211. Based on the interaction position of the γ photon passing through the aperture collimator with the crystal in the aperture coincidence event, calculate the coordinates of the center point of the corresponding crystal unit, and establish the straight line equation of the back projection line pointing from the center point of the crystal unit to the center point of the aperture collimator. S212. Determine the vector of the centerline of the slit collimator; S213. Based on the vector of the center line of the slit collimator, the coordinates of the crystal unit center, and the coordinates of the center point of the slit collimator, establish the plane equation of the back projection surface; S214. Based on the plane equation of the back-projection surface and the straight line equation of the back-projection line, solve for the three-dimensional coordinates of their intersection point. .

5. The method for improving the reconstruction speed and signal-to-noise ratio of a cascaded gamma imaging system according to claim 4, characterized in that, The step of mapping the intersection coordinates of the aperture coincidence event pairs in the field of view space to the corresponding field of view voxel numbers includes: S221. Pre-define the field of view to be reconstructed and divide it into a regular grid in three-dimensional space. Each grid cell is called a voxel. Let the FOV be... X The boundary of the direction is The number of voxels is ; Y Directional boundary is The number of voxels is The boundary in the Z direction is The number of voxels is The dimensions of each voxel in the three dimensions are: ； In the formula, , , They are respectively the voxels along the field of view X , Y , Z Dimensions in three directions; S222, Calculate the intersection point The voxel number corresponding to the intersection point is obtained by dividing the offset relative to the lower boundary of the FOV by the voxel size and rounding down. ； In the formula, This represents the floor function. , , They are respectively X , Y , Z The voxel index in the direction, and satisfying , , ; S223. For intersections that fall outside the FOV boundary, the corresponding events are removed.

6. The method for improving the reconstruction speed and signal-to-noise ratio of a cascaded gamma imaging system according to claim 1, characterized in that, The process of defining the region of interest centered on the voxel position represented by the voxel number and determining the boundary of the region of interest includes: S31. Taking the voxel position represented by the voxel number as the center, define a preset size within its neighborhood. m × n × l The voxel region is taken as the region of interest, where, m , n , l respectively along X axis, Y axis, Z Number of voxels along the axial direction; S32. Determine the boundary of the region of interest in the overall field of view voxel space based on the voxel space of the overall imaging field of view, the voxel number corresponding to the intersection point, and the size of the region of interest.

7. The method for improving the reconstruction speed and signal-to-noise ratio of a cascaded gamma imaging system according to claim 6, characterized in that, The step of determining the boundary of the region of interest in the overall field of view voxel space based on the voxel space of the overall imaging field of view, the voxel number corresponding to the intersection point, and the size of the region of interest includes: Assume the entire imaging field of view voxel space is in X axis, Y axis, Z The voxel index ranges along the axis are as follows: X axis: ,common Individual factors; Y axis: ,common Individual factors; Z axis: ,common Individual factors; Calculate the number of voxels extending outwards from the central voxel index in each dimension, and define the half-width in each direction as: ； In the formula, , , The regions of interest are respectively in X , Y , Z Half width in the direction; The theoretical voxel boundary of the region of interest is then: ； ； In the formula, , , The regions of interest are respectively in X , Y , Z Theoretical starting index in the direction; , , The regions of interest are respectively in X , Y , Z Theoretical termination index in direction, for Central voxel numbering; The obtained theoretical voxel boundaries are corrected; the correction formula is as follows: ； ； After correction, the region of interest is in X Axis effective voxel index range: Region of interest is Y Axis effective voxel index range: Region of interest is Z Axis effective voxel index range: .

8. The method for improving the reconstruction speed and signal-to-noise ratio of a cascaded gamma imaging system according to claim 1, characterized in that, Image iterative reconstruction is performed using the ordered subset expectation-maximization algorithm, with the following formula: ； In the formula, Voxel values ​​for the field of view and Global numbering of field voxels. For the number of iterations, It is the event number. The total number of events, The total number of voxels in the field of view. For the first The cascaded gamma photons produced by nuclear decay in the voxel of the field of view are the first The probability of an event. For the first The cascaded gamma photons produced by nuclear decay in the voxel of the field of view are the first The probability of an event. For the first A subset of data, superscript For the first In the first iteration Data subsets Iteration, It is the projected value. , , Indicates the region of interest.

9. A system for improving the reconstruction speed and signal-to-noise ratio of a cascaded gamma imaging system, wherein the system is used to implement the method for improving the reconstruction speed and signal-to-noise ratio of a cascaded gamma imaging system as described in any one of claims 1-8, characterized in that, The system includes: The subset partitioning module is used to encode the detector module based on the collected coincidence event data, combined with the gamma photon interaction information and collimator type, and to divide the coincidence event data into multiple subsets based on the energy information and the encoding. The voxel mapping module is used to perform direct back projection on the coincidence events in each subset to determine the projection intersection points and map them to the corresponding voxel numbers; The region determination module is used to delineate the region of interest centered on the voxel position represented by the voxel number, and to determine the boundary of the region of interest; The iterative reconstruction module is used to calculate the system transfer matrix corresponding to each pair of aperture coincidence events in each data subset. Based on the partitioned data subsets and the system transfer matrix, the ordered subset expectation maximization algorithm is used for iterative image reconstruction. When calculating the system transfer matrix, only the response probability between the coincidence event and each voxel in the corresponding region of interest is calculated. During the iterative reconstruction process, when calculating the projection and back projection values ​​of each pair of coincidence events, only the voxels in the region of interest are traversed.

Images