A compton imaging iterative reconstruction method for fast localization of radiation hotspots

Abstract

The application discloses a kind of Compton imaging iterative reconstruction methods for radiation hotspot fast positioning.This method is:1) from the coincidence event of detection screening imaging event;2) calculate the Compton conical projection of imaging event in projection space;3) according to each Compton conical projection generates a response matrix;4) projection space is divided into low resolution imaging space;5) after dividing into multiple subsystem response matrix, response matrix;Low resolution system response matrix is generated, and the integral value of its element corresponds to the subsystem response matrix;6) the iterative reconstruction of low resolution imaging space is carried out to obtain low resolution image;7) according to the maximum value index of the pixel of low resolution image corresponding subsystem response matrix, local high resolution iterative reconstruction is carried out, and the real angular orientation of radiation hotspot is obtained;8) according to the pixel of high resolution iterative reconstruction result of heat value area in low resolution image, projection space pixel is sequentially valued, and complete high resolution image is obtained.

Description

Technical Field

[0001] This invention is applied to the field of radiation imaging of gamma (γ) ray sources, and optimizes and improves the imaging algorithm of the Compton camera. In particular, it relates to an iterative reconstruction method for Compton imaging for rapid localization of radiation hotspots. Background Technology

[0002] In the field of radiation detection applications, advanced radiation detection technologies can help people perceive radiation and avoid harm. Radiation imaging has now become the mainstream technology in the field of nuclear radiation detection. It can vividly indicate the spatial distribution of radiation hotspots, visualizing the distribution of radioactive materials as hotspots in a "photographic" manner, and can accurately indicate the location of hotspots by combining optical video scenes. Compton imaging utilizes the Compton scattering effect when gamma rays interact with matter to image and locate radioactive sources. It uses a position-sensitive detector with multi-pixel resolution to acquire Compton scattering coincidence events, and then performs image reconstruction.

[0003] Commonly used Compton imaging image reconstruction algorithms mainly include Direct Backprojection (SBP) and Maximum Likelihood Expectation Maximization (MLEM). SBP is fast but suffers from poor angular resolution, while MLEM offers high image quality and angular resolution. However, MLEM is computationally intensive and time-complex, making it difficult to perform real-time image reconstruction or rapid imaging. Significant improvements have been made to accelerate MLEM, such as GPU-accelerated parallel computing. However, this requires more expensive hardware support, and some portable, miniaturized imaging systems, such as tablets or mobile phones, lack GPUs or parallel computing frameworks. Summary of the Invention

[0004] Given that current Compton imaging iterative reconstruction algorithms are computationally time-consuming and struggle to quickly locate radiation hotspots, the purpose of this invention is to provide a Compton imaging iterative reconstruction method for rapid location of radiation hotspots. This invention improves the algorithm's execution speed by partitioning and reorganizing the system response matrix required by the algorithm; it also boasts higher computational efficiency and lower computational time complexity compared to existing iterative reconstruction algorithms, enabling rapid location of radiation hotspots.

[0005] The technical solution of this invention is as follows:

[0006] A Compton imaging iterative reconstruction method for rapid localization of radiation hotspots includes the following steps:

[0007] 1) Record the matching events detected by the Compton camera in chronological order, and select the imaging events from the matching events;

[0008] 2) Divide the projection space of the Compton conic projection to obtain a two-dimensional projection array of size P×Q; calculate the Compton conic projection of the corresponding imaging instance in the projection space based on the position and energy information of the scattering and absorption points in each imaging instance; where P is the total number of longitudinal divisions of the projection space and Q is the total number of transverse divisions of the projection space.

[0009] 3) Generate an original high-resolution system response matrix based on the Compton conic projection corresponding to each imaging instance;

[0010] 4) Divide the projection space into a low-resolution imaging space represented by a two-dimensional array of size M×N; where M < P,

[0011] N < Q;

[0012] 5) Divide the original high-resolution system response matrix into M×N subsystem response matrices of equal size, where each subsystem response matrix has a size of (P / M)×(Q / N); generate a low-resolution system response matrix t of size M×N based on the integral values ​​of each subsystem response matrix. is The low-resolution system response matrix t is Each element in the matrix corresponds to the integral value of the response matrix of the subsystem.

[0013] 6) The low-resolution system response matrix t is As parameters input into the maximum likelihood expectation maximization iterative algorithm, a low-resolution image O is obtained through iterative reconstruction. low ;

[0014] 7) Based on the low-resolution image O low Find the corresponding subsystem response matrix A at the pixel position corresponding to the maximum pixel value in the matrix. k The subsystem response matrix A k Substituting this into the maximum likelihood expectation maximization iterative algorithm, an iterative reconstruction yields a local high-resolution image O. res Based on the local high-resolution image O res The true angular orientation of the radiating hotspot is obtained by calculating the maximum pixel value.

[0015] 8) From the low-resolution image O low Select a calorific value region and iteratively reconstruct it;

[0016] 9) Assign values ​​to the corresponding pixels in the projection space based on the pixel values ​​of each pixel in the reconstruction result of step 8) to obtain a complete high-resolution image.

[0017] Furthermore, the heat value region is the low-resolution image O. lowThe pixel values ​​in the image are greater than those in the low-resolution image O. low The area where the maximum pixel value is 50%.

[0018] Furthermore, in step 9), for the low-resolution image O low For pixels that have not undergone high-resolution reconstruction, their values ​​are divided by r×r to obtain the corresponding pixel values ​​in the projection space, and the complete high-resolution image is updated; where r is the set angular resolution of the low-resolution imaging space.

[0019] Furthermore, the projection space is divided into 180 degrees vertically and 360 degrees horizontally, with each degree as an interval, to obtain a two-dimensional projection array of size 180×360.

[0020] Furthermore, the true angular orientation of the radiation hotspot. Where θ is the longitudinal azimuth angle of the projection space. O is the lateral azimuth angle of the projection space. low_max .x represents a low-resolution image O low The x-coordinate of the maximum pixel value, O low_max .y represents a low-resolution image O low The ordinate of the maximum pixel value, O res_max .x represents a local high-resolution image O res The x-coordinate of the maximum pixel value, O res_max .y represents a local high-resolution image O res The ordinate of the maximum pixel value, r is the set angular resolution of the low-resolution imaging space.

[0021] The flowchart of the Compton imaging iterative reconstruction method for rapid localization of radiation hotspots in this invention is as follows: Figure 1 As shown, the specific implementation steps include:

[0022] 1. Record the coincidence events detected by the Compton camera in chronological order, select the imaging conditions for each coincidence event using the energy selection method, and count all imaging events.

[0023] 2. Based on the positional energy information of the scattering and absorbing points in each imaging event, calculate the Compton conic projection for the corresponding imaging event. The projection space is a 4π spherical field of view, and the projection space is determined according to the longitudinal azimuth angle θ and the lateral azimuth angle θ. It can be divided into 180 degrees vertically and 360 degrees horizontally. When calculating the projection array, the probability value under each azimuth angle can be calculated according to the array with each degree interval being 180×360.

[0024] 3. Calculate the Compton conic projection corresponding to each imaging instance. The set of all Compton conic projections is the original high-resolution system response matrix.

[0025] 4. Divide the original high-resolution imaging space (i.e., projection space) with an array size of 180×360 into a low-resolution imaging space of M×N, where M<180 and N<360.

[0026] 5. Following the imaging space partitioning method in step 4, the original high-resolution system response matrix is ​​divided into M×N subsystem response matrices, each with the same size (180 / M)×(360 / N). A new low-resolution system response matrix of size M×N is defined, with its elements corresponding to the integral values ​​of each subsystem response matrix.

[0027] 6. Input the low-resolution system response matrix as a parameter into the maximum likelihood expectation maximization iterative algorithm (MLEM) to perform iterative reconstruction of the low-resolution imaging space. The number of iterations is n, and the low-resolution image is obtained.

[0028] 7. Find the corresponding subsystem response matrix based on the position of the pixel with the maximum pixel value, substitute the subsystem response matrix into the MLEM algorithm for iterative calculation, and the number of iterations is also n to obtain a local high-resolution image; calculate the true angular orientation of the radiation hotspot based on the maximum pixel value in the reconstruction result (i.e., the local high-resolution image).

[0029] 8. Perform iterative reconstruction on pixels (hot regions) in low-resolution images whose pixel values ​​are greater than 50% of the maximum value.

[0030] 9. Based on the reconstruction results of each pixel in the low-resolution image, assign values ​​to the original high-resolution imaging spatial pixels of 180×360 sequentially to obtain the complete high-resolution image.

[0031] Compared with the prior art, the positive effects of the present invention are as follows:

[0032] 1) The proposed Compton imaging iterative reconstruction method for rapid localization of radiation hotspots improves the execution speed of the entire iterative algorithm while maintaining the high angular resolution of the original iterative algorithm, thus achieving rapid localization of radiation hotspots.

[0033] 2) The Compton imaging iterative reconstruction method for rapid localization of radiation hotspots proposed in this invention is an improvement on the original mature MLEM iterative reconstruction algorithm. It does not introduce additional variables, but only performs partitioning and recombination operations on the system response matrix in the prior parameters. The process is clear and concise.

[0034] As a supplement to the Compton imaging method, this invention proposes a fast and efficient iterative reconstruction method based on the partitioning and recombination of the system response matrix. Compared with the original Compton imaging iterative reconstruction algorithm, this method requires less computational time complexity and can achieve rapid localization of radiation hotspots. Attached Figure Description

[0035] Figure 1 This is a flowchart of the method of the present invention.

[0036] Figure 2 This is a conceptual diagram of the detector, scattering point, and absorption point in Compton imaging.

[0037] Figure 3 This is a calculation principle and projection diagram of the Compton cone.

[0038] Figure 4 A schematic diagram for generating the original high-resolution system response matrix.

[0039] Figure 5 This is a schematic diagram of the original high-resolution imaging space division.

[0040] Figure 6 Generate a schematic diagram for dividing the subsystem response matrix.

[0041] Figure 7 This is a schematic diagram illustrating the generation principle of the response matrix for a low-resolution system.

[0042] Figure 8 This is a schematic diagram of a low-resolution image.

[0043] Figure 9 This is a schematic diagram of the subsystem response matrix indexed based on the maximum pixel value in a low-resolution image.

[0044] Figure 10 This is a schematic diagram of the local high-resolution iterative reconstruction results.

[0045] Figure 11 This is a schematic diagram of the complete high-resolution image. Detailed Implementation

[0046] The present invention will now be described in further detail with reference to the accompanying drawings. The examples given are only for explaining the present invention and are not intended to limit the scope of the present invention.

[0047] The process of this invention is as follows Figure 1 As shown, the steps include:

[0048] 1. Record the coincidence events detected by the Compton camera in chronological order, and perform imaging condition selection for each coincidence event. Statistically analyze all imaging events. The concepts and principles of detectors, scattering points, and absorption points in Compton imaging are as follows: Figure 2 As shown.

[0049] 1.1 Select a typical Compton camera, or use a physical simulation tool to model and simulate the physical detection process of gamma rays based on the detector structure characteristics of the physical camera.

[0050] 1.2 When the same gamma ray hits two pixels of the detector, it is considered a coincidence case, and the position and deposited energy of each pixel are recorded.

[0051] 1.3 Based on the energy values ​​of the two pixels, the imaging conditions are determined by the energy selection method to distinguish the scattering point and absorption point corresponding to the two pixels respectively.

[0052] 2. Based on the positional energy information of the scattering and absorbing points in each imaging event, calculate its Compton conic projection, with a 4π spherical field of view in the projection space, according to the azimuth angle θ and It can be divided into 180 degrees vertically and 360 degrees horizontally. When calculating the projection array, the probability value for each azimuth angle can be calculated according to an array with each degree interval being 180×360, such as... Figure 3 As shown.

[0053] 2.1 The position and energy information of the scattering point and absorption point of the i-th Compton coincidence imaging case are (x1, y1, z1, E1) and (x2, y2, z2, E2), respectively, where (x1, y1, z1) is the position information of the scattering point, E1 is the energy information of the scattering point, p1 is the position of the scattering point (x1, y1, z1), p2 is the position of the absorption point (x2, y2, z2), and n i Let θ be the unit vector of the Compton cone axis, and let θ be the Compton scattering angle. i The expression is as follows:

[0054]

[0055] Where, subscript i represents the i-th Compton coincidence imaging instance, m e c 2 This is the static electron energy.

[0056] 2.2 The calculation principle and projection of the Compton cone are as follows: Figure 3 As shown, if λ j Let represent the probability value of gamma rays emitted from voxel j on the imaging sphere. This is calculated iteratively for each voxel in the 4π imaging space, yielding a Compton's cone for each imaging scattering event, represented on the sphere as a ring of probability values. The 4π spherical field of view is calculated according to the azimuth angle θ and... It can be divided into 180 degrees vertically and 360 degrees horizontally. When calculating the projection array, the probability value for each azimuth angle can be calculated using an array with each degree interval being 180×360. Assume the angle error is Δθ. i The value is typically taken as 10°. The angle between the orientation vector of voxel j and the unit vector ni of the reconstructed cone axis is denoted as α. ij The probability value of the Compton conical projection is expressed as:

[0057]

[0058] 3. Calculate the Compton conic projection for each event; the set of all projections is the original high-resolution system response matrix T. ij Let 'i' represent the i-th imaging event, and 'j' represent the j-th pixel in the 180°×360° imaging space. The number of 'j' pixels ranges from 1 to 180×360 (i.e., pixels are divided into 1° intervals). Assuming there are a total of I imaging events, the matrix size is I×180×360. The principle of the original high-resolution system response matrix is ​​as follows: Figure 4 As shown.

[0059] 4. Divide the original high-resolution imaging space of size 180×360 into a low-resolution imaging space of size M×N, where M < 180 and N < 360. Assuming M is 6 and N is 12, this represents an angular resolution r = 30°. Figure 5 As shown.

[0060] 5. Following the imaging space partitioning method in step 4, the original high-resolution system response matrix is ​​divided into M×N (6×12) subsystem response matrices A. k The range of k is 1 to 6 × 12. The principle for generating the subsystem response matrix is ​​as follows: Figure 6 As shown, the array of each original Compton cone is divided into M×N (6×12) parts. The parts with consistent coordinates are then combined to form the subsystem response matrix. There are M×N (6×12) subsystem response matrices, each with a size of I×30×30. A new low-resolution system response matrix t is defined. is Where the subscript i represents the i-th Compton cone, s represents the coordinates of the corresponding pixel in the low-resolution imaging space, and its value is the integral value of the corresponding subsystem response matrix (i.e., the sum of the matrix array). The generation principle of the low-resolution system response matrix is ​​as follows: Figure 7 As shown. The response matrix A of the first subsystem... k＝1 For example, it is equivalent to θ = 0° to 30° in the original high-resolution system response matrix and A small portion of it, its matrix size is I×30×30, and the corresponding low-resolution system response matrix t is＝1 For A k＝1 The sum of the values ​​of A k＝1 The mathematical form is:

[0061]

[0062] 6. The low-resolution system response matrix t is Substituting the values ​​into the MLEM algorithm, iterative reconstruction of the low-resolution imaging space is performed, with n = 10 iterations, to obtain the low-resolution image, such as... Figure 8As shown, the time complexity of the algorithm is O(M×N=6×12), and the iteration formula is:

[0063]

[0064] 7. Acquire low-resolution images O low The maximum value of the pixels in the range, and based on its pixel position O low_max The subsystem response matrix A corresponding to the index (x,y)=(2,5) is A k＝12×2+5＝29 ,like Figure 9 As shown, local high-resolution iterative reconstruction is performed. The iterative formula is the same as the above formula, only needing to be replaced with the subsystem response matrix A. k＝29 The number of iterations is n=10, and the reconstruction result is O. res like Figure 10 As shown, the time complexity of the algorithm is O(r×r=30×30). The true angular orientation of the radiation hotspot is calculated based on the maximum pixel value in the reconstruction result. Considering the algorithm execution time in step 6, the total time complexity of the algorithm is O(M×N+r×r=972). The time complexity required for the original iterative algorithm is O(180×360=64800), which is 64800 / 972=66 times faster.

[0065] Figure 10 The coordinate of the maximum value in is O res_max (x,y)=(5,16), when performing low-resolution processing of the system response matrix, the angular resolution r=30°, therefore, the true angular orientation of the radiation hotspot is...

[0066] 8. Based on the reconstruction results of the low-resolution imaging space O low ( Figure 8 As shown in the figure, pixels with values ​​greater than 50% of the peak value are selected as hotspot region pixels, and local high-resolution iterative reconstruction is performed on them. Then, the original 180×360 high-resolution imaging spatial pixels are assigned values ​​sequentially. For low-resolution image O low The non-thermal regions in the image are divided by 900 (30×30) to obtain the corresponding pixel values ​​in the original high-resolution imaging space, ultimately yielding the complete high-resolution image, as shown below. Figure 11 As shown.

[0067] In summary, the above are merely preferred embodiments of the present invention and are not intended to limit the scope of protection of the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A Compton imaging iterative reconstruction method for rapid localization of radiation hotspots, comprising the following steps: 1) Record the matching events detected by the Compton camera in chronological order, and select the imaging events from the matching events; 2) Divide the projection space of the Compton conic projection to obtain a two-dimensional projection array of size P×Q; calculate the Compton conic projection of the corresponding imaging instance in the projection space based on the position and energy information of the scattering and absorption points in each imaging instance; where P is the total number of longitudinal divisions of the projection space and Q is the total number of transverse divisions of the projection space. 3) Generate an original high-resolution system response matrix based on the Compton conic projection corresponding to each imaging instance; 4) Divide the projection space into a low-resolution imaging space represented by a two-dimensional array of size M×N; where M < P, N < Q; 5) Divide the original high-resolution system response matrix into M×N subsystem response matrices of equal size, where each subsystem response matrix has a size of (P / M)×(Q / N); generate a low-resolution system response matrix t of size M×N based on the integral values ​​of each subsystem response matrix. is The low-resolution system response matrix t is Each element in the matrix corresponds to the integral value of the response matrix of the subsystem. 6) The low-resolution system response matrix t is As parameters input into the maximum likelihood expectation maximization iterative algorithm, a low-resolution image O is obtained through iterative reconstruction. low ; 7) Based on the low-resolution image O low Find the corresponding subsystem response matrix A at the pixel position corresponding to the maximum pixel value in the matrix. k The subsystem response matrix A k Substituting this into the maximum likelihood expectation maximization iterative algorithm, an iterative reconstruction yields a local high-resolution image O. res Based on the local high-resolution image O res The true angular orientation of the radiating hotspot is obtained by calculating the maximum pixel value. 8) From the low-resolution image O low Select a calorific value region and iteratively reconstruct it; 9) Assign values ​​to the corresponding pixels in the projection space based on the pixel values ​​of each pixel in the reconstruction result of step 8) to obtain a complete high-resolution image.

2. The method according to claim 1, characterized in that, The heat value region is the low-resolution image O. low The pixel values ​​in the image are greater than those in the low-resolution image O. low The area where the maximum pixel value is 50%.

3. The method according to claim 1, characterized in that, In step 9), for the low-resolution image O low For pixels that have not undergone high-resolution reconstruction, their values ​​are divided by r×r to obtain the corresponding pixel values ​​in the projection space, and the complete high-resolution image is updated; where r is the set angular resolution of the low-resolution imaging space.

4. The method according to claim 1, characterized in that, The projection space is divided into 180 degrees vertically and 360 degrees horizontally, with each degree as an interval, to obtain a two-dimensional projection array of size 180×360.

5. The method according to claim 1, characterized in that, True angular location of radiation hotspots Where θ is the longitudinal azimuth angle of the projection space. O is the lateral azimuth angle of the projection space. low_max .x represents a low-resolution image O low The x-coordinate of the maximum pixel value, O low_max .y represents a low-resolution image O low The ordinate of the maximum pixel value, O res_max .x represents a local high-resolution image O res The x-coordinate of the maximum pixel value, O res_max .y represents a local high-resolution image O res The ordinate of the maximum pixel value, r is the set angular resolution of the low-resolution imaging space.

Images