Compton imaging iteration reconstruction method used for radiation hot spot rapid positioning

Abstract

The invention discloses a Compton imaging iteration reconstruction method used for radiation hot spot rapid positioning. The method comprises the following steps of: 1) screening an imaging case from detected coincidence cases; 2) calculating Compton cone projection of the imaging case in a projection space; 3) generating a response matrix according to each Compton cone projection; 4) dividing the projection space into a low-resolution imaging space; 5) dividing the response matrix into a plurality of subsystem response matrixes; generating a low-resolution system response matrix, wherein elements of the low-resolution system response matrix correspond to integral values of the subsystem response matrix; 6) performing iterative reconstruction of the low-resolution imaging space to obtain a low-resolution image; 7) according to the sub-system response matrix corresponding to the pixel maximum value index of the low-resolution image, performing local high-resolution iterative reconstruction to obtain the real angular orientation of the radiation hotspot; and 8) according to the high-resolution iteration reconstruction result of the pixels of the calorific value area in the low-resolution image, assigning values to the pixels of the projection space in sequence to obtain a complete high-resolution image.

Description

Technical Field

[0001] The present invention is applied to the field of radiation imaging of gamma (γ) ray sources, optimizes and improves the imaging algorithm of a Compton camera, and particularly relates to a Compton imaging iterative reconstruction method for rapid positioning of radiation hot spots. Background Art

[0002] In the field of radiation detection applications, advanced radiation detection technology 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, and use the "photographing" method to visualize the distribution of radioactive materials in the form of hotspots, and can accurately give the distribution location of hotspots in conjunction with optical video scenes. Compton imaging uses the Compton scattering effect when gamma rays interact with matter to image and locate the radiation source, and obtains Compton scattering coincidences through a position-sensitive detector with multi-pixel resolution, and then reconstructs the image.

[0003] The commonly used Compton imaging image reconstruction algorithms mainly include the direct back-projection reconstruction algorithm (SBP) and the maximum likelihood expectation maximization iterative reconstruction algorithm (MLEM). Among them, the direct back-projection reconstruction algorithm is fast, but the angular resolution of the reconstructed image is poor, while the MLEM iterative reconstruction algorithm has good imaging quality and high angular resolution. However, the MLEM iterative reconstruction algorithm has a large amount of calculation and high computational time complexity, making it difficult to perform real-time image reconstruction or fast imaging. For the accelerated execution of the MLEM iterative reconstruction algorithm, some studies have made significant improvements, such as parallel computing based on GPU acceleration, but the required hardware support cost is higher. In some portable and miniaturized imaging systems, there is no GPU or computing framework that can perform parallel computing, such as tablets or mobile phones. Summary of the invention

[0004] In view of the fact that the current iterative reconstruction algorithm of Compton imaging has a long calculation time and is difficult to quickly locate radiation hot spots, the purpose of the present invention is to provide a Compton imaging iterative reconstruction method for quickly locating radiation hot spots. The present invention improves the execution speed of the algorithm by dividing and reorganizing the system response matrix required by the algorithm; the present invention has higher calculation efficiency and lower calculation time complexity than the original iterative reconstruction algorithm, and can realize the rapid location of radiation hot spots.

[0005] The technical solution of the present invention is:

[0006] A Compton imaging iterative reconstruction method for rapid positioning of radiation hot spots, comprising the following steps:

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

[0008] 2) Dividing the projection space of the Compton cone projection to obtain a two-dimensional projection array of size P×Q; calculating the Compton cone projection of the corresponding imaging case in the projection space according to the position energy information of the scattering point and the absorption point in each imaging case; wherein P is the total number of longitudinally divided spaces of the projection space, and Q is the total number of horizontally divided spaces of the projection space;

[0009] 3) generating an original high-resolution system response matrix according to the Compton conic projection corresponding to each imaging case;

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

[0011] N<Q;

[0012] 5) Dividing the original high-resolution system response matrix into M×N subsystem response matrices of the same size, the size of the subsystem response matrix is ​​(P / M)×(Q / N); generating a low-resolution system response matrix t of size M×N according to the integral value of each subsystem response matrix is , where the low-resolution system response matrix t is Each element in corresponds to an integral value of the subsystem response matrix;

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

[0014] 7) According to the low-resolution image O low The corresponding subsystem response matrix A is found by taking the pixel position corresponding to the maximum pixel value in k , the subsystem response matrix A k Substitute it into the maximum likelihood expectation maximization iterative algorithm and iteratively reconstruct a local high-resolution image O res ; Based on the local high-resolution image O res The true angular position of the radiation hotspot is calculated by using the maximum value of the pixels in;

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

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

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

[0018] Further, in step 9), for the low-resolution image O low For pixels that have not been reconstructed with high resolution, their values ​​are divided by r×r as the pixel values ​​of the corresponding pixels in the projection space to update the complete high-resolution image; wherein r is the angular resolution of the set low-resolution imaging space.

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

[0020] Furthermore, the true angular position of the radiation hotspot Among them, θ is the longitudinal azimuth angle of the projection space, is the transverse azimuth of the projection space, O low_max .x is a low-resolution image O low The horizontal coordinate of the maximum value of the pixel in O low_max .y is the low-resolution image O low The vertical coordinate of the maximum value of the pixel in O res_max .x is the local high-resolution image O res The horizontal coordinate of the maximum value of the pixel in O res_max .y is the local high-resolution image O res The vertical coordinate of the maximum value of the pixel in the image is r, and r is the angular resolution of the set low-resolution imaging space.

[0021] The flowchart of the Compton imaging iterative reconstruction method for rapid positioning of radiation hot spots of the present invention is as follows: Figure 1 As shown, the specific implementation steps include:

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

[0023] 2. According to the position energy information of the scattering point and the absorption point in each imaging case, the Compton cone projection of the corresponding imaging case is calculated. The projection space is a 4π spherical field of view. The projection space is divided into the longitudinal azimuth angle θ and the transverse azimuth angle It can be divided into 180 degrees vertically and 360 degrees horizontally. When calculating the projection array, the probability value at each azimuth angle can be calculated according to the array with each degree interval of 180×360.

[0024] 3. Count the Compton conic projections corresponding to each imaging case. 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. According to the imaging space division method in step 4, the original high-resolution system response matrix is ​​divided into M×N subsystem response matrices, and the size of each subsystem response matrix is ​​the same, which is (180 / M)×(360 / N). Define a new low-resolution system response matrix with a size of M×N, whose elements correspond to the integral value of each subsystem response matrix.

[0027] 6. The low-resolution system response matrix is ​​input 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 a low-resolution image is obtained.

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

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

[0030] 9. According to the reconstruction results of each pixel in the low-resolution image, the original high-resolution imaging space pixels of 180×360 are assigned values ​​in turn to obtain a complete high-resolution image.

[0031] Compared with the prior art, the present invention has the following positive effects:

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

[0033] 2) The Compton imaging iterative reconstruction method for rapid positioning of radiation hotspots proposed in the present invention is an improvement on the original mature MLEM iterative reconstruction algorithm. It does not introduce additional variables and only divides and reorganizes the system response matrix in the prior parameters. The process is clear and concise.

[0034] As a supplement to the Compton imaging method, the present invention proposes a fast and efficient iterative reconstruction method based on the partitioning and reorganization of the system response matrix. Compared with the original Compton imaging iterative reconstruction algorithm, this method requires lower computational time complexity and can achieve rapid positioning of radiation hotspots. BRIEF DESCRIPTION OF THE DRAWINGS

[0035] Figure 1 The figure is a flow chart of the method of the present invention.

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

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

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

[0039] Figure 5 Schematic diagram of the original high-resolution imaging space division.

[0040] Figure 6 Generate schematics for partitioning of subsystem response matrices.

[0041] Figure 7 Schematic diagram for generating the response matrix of a low-resolution system.

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

[0043] Fig. 9 Schematic diagram of the subsystem response matrix corresponding to the maximum pixel value index in the low-resolution image.

[0044] Fig.10 Schematic diagram of local high-resolution iterative reconstruction results.

[0045] Fig.11 A schematic diagram of a complete high-resolution image. DETAILED DESCRIPTION

[0046] The present invention is further described in detail below in conjunction with the accompanying drawings. The examples given are only used to explain the present invention but not to limit the scope of the present invention.

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

[0048] 1. Record the matching events detected by the Compton camera in chronological order, and select the imaging conditions for each matching event, and count all imaging events. The concept and principle of detectors, scattering points, and absorption points in Compton imaging are as follows: Figure 2 shown.

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

[0050] 1.2 When the same gamma ray hits two pixels of the detector, it is considered a coincidence event 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 judged by the energy selection method to distinguish the scattering points and absorption points corresponding to the two pixels.

[0052] 2. According to the position energy information of the scattering point and the absorption point in each imaging case, the Compton conic projection is calculated. The projection space is a 4π spherical field of view. It can be divided into 180 degrees vertically and 360 degrees horizontally. When calculating the projection array, the probability value at each azimuth angle can be calculated according to the array with an interval of 180×360 for each degree, such as Figure 3 shown.

[0053] 2.1 The position and energy information of the scattering point and absorption point of the i-th Compton coincidence imaging event 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), n i is the unit vector of the Compton cone axis, and its Compton scattering angle θ i It is expressed as:

[0054]

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

[0056] 2.2 Calculation principle and projection of Compton cone Figure 3 As shown, if λ j Represents the probability value of gamma rays emitted from voxel j on the imaging sphere. Each voxel in the 4π imaging space is traversed and calculated to obtain the Compton cone of each imaging scattering event, which is represented as a probability value ring on the sphere. The 4π spherical field of view, 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 at each azimuth angle can be calculated according to the array with an interval of 180×360 for each degree. Assume that the angle error is Δθ i , which is generally set to 10°. The angle between the orientation vector of voxel j and the unit vector ni of the reconstruction cone axis is denoted as α ij , then the probability value of the Compton conic projection is expressed as:

[0057]

[0058] 3. Count the Compton cone projections corresponding to each case, and the collection of all projections is the original high-resolution system response matrix T ij , i represents the i-th imaging case, j represents the j-th pixel in the 180°×360° imaging space, the number of j is 1 to 180×360 (i.e., pixels are divided into 1° intervals), assuming that there are a total of I imaging cases, then the matrix size is I×180×360, the principle of the original high-resolution system response matrix is ​​as follows Figure 4 shown.

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

[0060] 5. According to the imaging space division method in step 4, the original high-resolution system response matrix is ​​divided into M×N (6×12) subsystem response matrices A k , k ranges from 1 to 6×12, and the principle of partitioning and generating the subsystem response matrix is ​​as follows Figure 6 As shown, each original Compton cone array is divided into M×N (6×12) parts, and then the parts with consistent coordinates are combined to form the subsystem response matrix. There are M×N (6×12) subsystem response matrices, and the size of each subsystem response matrix is ​​I×30×30. Define the new low-resolution system response matrix t 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. Take the first subsystem response matrix A k＝1 For example, it is equivalent to the original high-resolution system response matrix θ = 0° ~ 30° and A small part of it, its matrix size is I×30×30, the corresponding low-resolution system response matrix t is＝1 A k＝1 The value and A k＝1 The mathematical form of is:

[0061]

[0062] 6. The low-resolution system response matrix t is Substitute it into the MLEM algorithm to iteratively reconstruct the low-resolution imaging space. The number of iterations is n = 10 times, and a low-resolution image is obtained, 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. Get low-resolution image O low The maximum value of the pixel in and according to its pixel position O low_max The subsystem response matrix A corresponding to the index (x,y) = (2,5) k＝12×2+5＝29 ,like Fig. 9 As shown, local high-resolution iterative reconstruction is performed. The iterative formula is the same as the above formula, and only needs to be replaced by the subsystem response matrix A k＝29 , the number of iterations is n = 10 times, and the reconstruction result is O res like Fig.10 As shown, the time complexity of the algorithm execution is O(r×r=30×30). The true angular orientation of the radiation hotspot is calculated according to the maximum pixel value in the reconstruction result. Taking into account the algorithm execution time of step 6, the total time complexity of the algorithm execution is O(M×N+r×r=972). The time complexity required for the original iterative algorithm execution is O(180×360=64800), which is 64800 / 972=66 times faster.

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

[0066] 8. Based on the reconstruction results of low-resolution imaging space O low ( Figure 8 As shown in Figure 2, pixels with values ​​greater than 50% of the peak are selected as hotspot pixels, and local high-resolution iterative reconstruction is performed on them respectively. Then, the original high-resolution imaging space pixels of 180×360 are assigned values ​​one by one. low The non-thermal value area in , its value is divided by 900 (30 × 30) as the pixel value of the corresponding area in the original high-resolution imaging space, and finally the complete high-resolution image is obtained, such as Fig.11 shown.

[0067] In summary, the above is only a preferred implementation method of the present invention, and is not intended to limit the protection scope of the present invention. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present invention shall be included in the protection scope of the present invention.

Claims

1. A Compton imaging iterative reconstruction method for rapid positioning of radiation hot spots, comprising the following steps: 1) Record the matching events detected by the Compton camera in chronological order, and select imaging events from the matching events; 2) Dividing the projection space of the Compton cone projection to obtain a two-dimensional projection array of size P×Q; calculating the Compton cone projection of the corresponding imaging case in the projection space according to the position energy information of the scattering point and the absorption point in each imaging case; wherein P is the total number of longitudinally divided spaces of the projection space, and Q is the total number of horizontally divided spaces of the projection space; 3) generating an original high-resolution system response matrix according to the Compton conic projection corresponding to each imaging case; 4) Dividing the projection space into a low-resolution imaging space represented by a two-dimensional array of size M×N; wherein M<P, N<Q; 5) Dividing the original high-resolution system response matrix into M×N subsystem response matrices of the same size, the size of the subsystem response matrix is ​​(P / M)×(Q / N); generating a low-resolution system response matrix t of size M×N according to the integral value of each subsystem response matrix is , where the low-resolution system response matrix t is Each element in corresponds to an integral value of the subsystem response matrix; 6) The low-resolution system response matrix t is As a parameter input into the maximum likelihood expectation maximization iterative algorithm, a low-resolution image O is obtained by iterative reconstruction. low ; 7) According to the low-resolution image O low The corresponding subsystem response matrix A is found by taking the pixel position corresponding to the maximum pixel value in k , the subsystem response matrix A k Substitute it into the maximum likelihood expectation maximization iterative algorithm and iteratively reconstruct a local high-resolution image O res ; Based on the local high-resolution image O res The true angular position of the radiation hotspot is calculated by using the maximum value of the pixels in; 8) From the low-resolution image O low Select a heat value area in and iteratively reconstruct it; 9) Assigning values ​​to corresponding pixels in the projection space according to the pixel value 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 area is the low-resolution image O low The pixel value in is greater than that 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 been reconstructed with high resolution, their values ​​are divided by r×r as the pixel values ​​of the corresponding pixels in the projection space to update the complete high-resolution image; wherein r is the angular resolution of the set low-resolution imaging space.

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

5. The method according to claim 1, characterized in that The true angular position of the radiation hotspot Among them, θ is the longitudinal azimuth angle of the projection space, is the transverse azimuth of the projection space, O low_max .x is a low-resolution image O low The horizontal coordinate of the maximum value of the pixel in O low_max .y is the low-resolution image O low The vertical coordinate of the maximum value of the pixel in O res_max .x is the local high-resolution image O res The horizontal coordinate of the maximum value of the pixel in O res_max .y is the local high-resolution image O res The vertical coordinate of the maximum value of the pixel in the image is r, and r is the angular resolution of the set low-resolution imaging space.

Images