Fast 3D detection imaging method for industrial gamma photon detector

Abstract

The application discloses a kind of fast 3D detection imaging methods for the region of interest of industrial gamma photon detector, the method comprises the following steps, according to the relative position of measured object and detection system, the profile position of measured object is calculated, the pixel corresponding in profile region is regarded as the region of interest. The probability contribution weight of each response line to the pixel in the region of interest is calculated by solid angle-distance model, and the system matrix of the region of interest is constructed. The system matrix of interest is calculated in parallel on GPU by CUDA architecture, and the image of the region of interest is iteratively reconstructed quickly. And, in each image iterative reconstruction, a Gaussian filter is added to further suppress the interference of internal cluttered scatter noise of image. The application calculates the system matrix for the region of interest, not only improves the imaging speed, and the solid angle-distance model makes the calculation of system matrix more accurate.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of image reconstruction and image processing technology, and in particular to a fast 3D detection imaging method for a region of interest of an industrial gamma photon detector. BACKGROUND

[0002] Positron Emission Tomography (PET) utilizes a pair of gamma photons generated by positron annihilation, the direction of the pair of gamma photons differs by 180°, has strong penetration and electrical neutrality, and is less affected by materials, structures and surrounding electric and magnetic fields, so the technology is widely used in the medical field. At present, the positron annihilation technology is also widely concerned in the field of industrial nondestructive testing. The high-energy gamma photon pair generated by positron annihilation can well penetrate composite materials and metal materials, and the positron nuclide is labeled on the working medium or marker and placed inside the object, so that the state information inside the object can be presented in the form of a three-dimensional image. As the gamma photon imaging technology tends to be mature, in order to adapt to more demand scenarios, the detector space is also getting larger and larger, but the slow imaging speed and the large amount of storage resources restrict the further development of positron technology, so it is of great significance to study the fast reconstruction of gamma photon images

[0003] List-mode data is a mode of recording the corresponding detection crystal number, the time of arrival at both ends of the crystal and the energy information of each coincidence response line (LOR) in a list, which can realize real-time and dynamic reconstruction; however, the current method of image reconstruction using list-mode data has a huge amount of calculation, which needs to calculate the probability weight of all pixels in the whole space corresponding to each LOR, so as to update the image, and how to ensure the image quality while quickly calculating or reducing the amount of calculation has always been a difficulty in gamma photon image reconstruction.

[0004] In addition, the traditional statistical iterative reconstruction algorithm considers that the image reconstruction is based on the assumption that the annihilation point is equally probable on the LOR, so all pixels within a certain length from the LOR are assigned a probability weight, which will cause the image to appear in places where there should be no radioactive material, the intensity of the places where the image should appear is weakened, and the pixel value of the image deviates from the actual value, not only reducing the radioactivity concentration in the active area, but also making the image size larger and the boundary blurred, and appearing the phenomenon of "expansion". How to enhance the imaging of the active area and reduce or even eliminate the imaging of the non-active area is also a difficulty in gamma photon image reconstruction. SUMMARY

[0005] The technical problems solved by the present application are to provide a fast 3D detection imaging method for an industrial gamma photon detector, and to accelerate the image reconstruction algorithm of list-mode data while improving the quality of the reconstructed image.

[0006] To solve the above technical problems, the technical solution adopted by the present application is a fast 3D detection imaging method for an industrial gamma photon detector, which specifically comprises the following steps:

[0007] Step 1: According to the relative position of the measured object and the detection system, the contour position of the measured object, i.e. the region of interest, is obtained, and the voxel information in the region of interest is recorded and stored;

[0008] The contour position of the measured object is obtained according to the position of the measured workpiece relative to the ring-shaped detector, the center of the measured workpiece is made to coincide with the center of the ring-shaped detector, the voxel position information in the contour activity region is calculated, and the specific operation is as follows: the size of the rectangle inscribed in the ring-shaped detector is A×B×C, the distance from the center of the measured workpiece to the contour edge is dx\dy\dz, the size of the reconstructed image is a×b×c, the distance from the center of the reconstructed image to the contour edge of the measured object to be reconstructed is tx / ty / tz, and the contour positions of the workpiece to be reconstructed along the x, y and z coordinate axes are calculated according to formulas (1), (2) and (3) to obtain the voxel position information in the contour activity region to be reconstructed:

[0009]

[0010] Step 2: The system matrix element values are divided into small units that can be executed independently, the calculation of the probability weight corresponding to each voxel in the region of interest by each response line is independent of each other, and the update of all voxels corresponding to each response line during image iterative reconstruction is also independent of each other, task allocation is performed through the CUDA architecture, each response line is allocated an ID of a block, the probability weight corresponding to each voxel in the region of interest is calculated in parallel, and the system matrix calculation in the region of interest is completed;

[0011] The calculation of the probability weight corresponding to each voxel in the region of interest is as follows: the included angle between the voxel center point and the detection crystal is a solid angle, the pair of detection crystals corresponding to the response line passing through the voxel are A and B, and the minimum solid angle value of the detection crystal pair to the voxel is taken as the weight of the system matrix in the image reconstruction, i.e. the element value a of the system matrix ij , and the probability weight of the voxel is μ×a ij , wherein the proportion coefficient μ of the weight is the distance ratio of the voxel to the detection crystal A and the detection crystal B, μ=D / P; D is the distance of the detection crystal A corresponding to the minimum solid angle to the voxel center, and P is the distance of the voxel center to the other detection crystal B.

[0012] Step 3: reassign a new ID of block to each response line, and complete 3D-OSEM image iterative reconstruction based on the region of interest system matrix, the iterative reconstruction formula is as follows:

[0013]

[0014] Wherein, S n is the nth subset, y i represents the ith LOR, λ j represents the jth voxel, and respectively represent the results of the reconstructed image in the iterative k+1 and k times, a ij is an element in the system matrix A, which represents the probability that a pair of gamma photons released by the jth voxel is captured by the crystal connected by the ith LOR. N represents the total number of voxels to be reconstructed, a il represents the probability that a pair of gamma photons released by the lth voxel is captured by the crystal connected by the ith LOR.

[0015] As a preferred scheme, in the image iterative reconstruction in step 3, a Gaussian filter is added in the iterative update.

[0016] As a preferred scheme, when the CUDA architecture is used to reconstruct the image, the resolution of the reconstructed image is MxM, MxM threads are opened in the block, the highest number of threads opened in each block is 1024, the lowest number of threads opened in each block is 32, and the number of threads opened must be a multiple of 32.

[0017] The beneficial effects of the present application are:

[0018] (1) According to the relative position of the measured object and the detection system, the position of the measured object contour edge can be accurately calculated, and the position in the contour is determined as the region of interest. Compared with the existing three-dimensional system matrix calculation method, the present application only needs to calculate the system matrix elements corresponding to the voxels in the region of interest, and does not need to calculate the voxels in a complete three-dimensional system matrix, which can greatly reduce the calculation amount and storage space of the three-dimensional system matrix in image reconstruction, and provides a feasible solution for fast image reconstruction of large space, large aperture and long axial industrial gamma photon detection system, and improves the production speed of industrial detection

[0019] (2) The application adds distance as a weight proportion coefficient on the basis of the traditional system matrix calculation model solid angle, and different voxels correspond to the system matrix according to the distance to perform different ratio distribution of the matrix weight value, so that the accuracy of the weight of a pair of gamma photons generated by annihilation of each voxel and detected by the response line can be improved. After the solid angle-distance model calculates the system matrix elements corresponding to the voxels in the region of interest, the accurate region of interest system matrix can be used to reconstruct a three-dimensional reconstruction image with clear edges and high resolution.

[0020] (3) The application adds a Gaussian filter in each image iterative reconstruction to suppress the internal cluttered scattered point noise interference of the image, so as to further improve the three-dimensional imaging quality and meet the requirements of industrial gamma photon detection systems in nondestructive testing.

[0021] The photon detection system is applied to the requirements of industrial nondestructive testing. BRIEF DESCRIPTION OF DRAWINGS

[0022] Figure 1 The flow chart of the region of interest fast 3D detection imaging method of the application;

[0023] Figure 2 The parallel calculation of the CUDA architecture on the GPU and the image iterative reconstruction flow chart;

[0024] Figure 3 The solid angle-distance model calculation schematic diagram;

[0025] Figure 4 The simulation model schematic diagram;

[0026] Figure 5 The traditional solid angle model, the solid angle-distance model of the application, and the image reconstruction result schematic diagram obtained after the Gaussian filter is added. DETAILED DESCRIPTION

[0027] The specific implementation of the application will be described in detail below with reference to the accompanying drawings.

[0028] As shown in the figure, the region of interest fast 3D detection imaging method for the industrial gamma photon detector specifically includes the following steps: Figures 1-3

[0029] Step 1: According to the relative position of the measured object and the detection system, the contour position of the measured object, i.e. the region of interest, is obtained, and the voxel information in the region of interest is recorded and stored;

[0030] ​The contour position of the object under test is obtained based on the position of the workpiece relative to the ring detector. The center of the workpiece is made to coincide with the center of the ring detector. The voxel position information in the contour activity area is calculated. Specifically, the rectangle inscribed in the ring detector is A×B×C, the distance from the center of the workpiece to the contour edge is dx\dy\dz, the size of the reconstructed image is a×b×c, and the distance from the center of the reconstructed image to the contour edge of the object to be reconstructed is tx / ty / tz. The contour position of the workpiece to be reconstructed along the coordinate axes x, y, z is calculated according to equations (1), (2), and (3) to obtain the voxel position information to be reconstructed in the contour activity area.

[0031]

[0032] Step 2: As Figure 2 As shown, the computation system matrix element values ​​are divided into small units that can be executed independently. The calculation of the probability weights corresponding to voxels in the region of interest for each type of response line is independent. During image iterative reconstruction, the updates of all voxels corresponding to each type of response line are also independent. Task allocation is performed through the CUDA architecture, assigning a block ID to each type of response line, and calculating the probability weights corresponding to voxels in the region of interest in parallel until the computation of the region of interest system matrix is ​​completed.

[0033] The CUDA architecture records the resolution of the reconstructed image as M×M, and opens M×M threads in the block. Each block opens a maximum of 1024 threads and a minimum of 32 threads, and the number of threads must be a multiple of 32.

[0034] like Figure 3 As shown, the probability weights corresponding to voxels within the region of interest are calculated as follows: the angle between the voxel center point and the detector crystal is the solid angle; the pair of detector crystals corresponding to the response lines passing through the voxel are A and B, respectively; and the minimum solid angle subtended by the detector crystals with respect to the voxel is calculated. Figure 3 In The angle value is used as the weight of the system matrix in image reconstruction, i.e., the element value a of the system matrix. ij Then the probability weight of the voxel is μ×a ij In the formula, the weighting coefficient μ is the ratio of the distance between the voxel and detector crystals A and B, μ=D / P; D is the distance between detector crystal A and the center of the voxel corresponding to the smallest solid angle, and P is the distance between the center of the voxel and another detector crystal B;

[0035] Step 3: reassign a new ID of block to each response line, and complete 3D-OSEM image iterative reconstruction based on the region of interest system matrix, and add a Gaussian filter in iterative update. The iterative reconstruction formula is as follows:

[0036]

[0037] wherein S n is the nth subset, y i represents the ith LOR, λ j represents the jth voxel, and respectively represent the results of the reconstructed image in the iterative k+1 and k times, a ij is an element in the system matrix A, and represents the probability that a pair of γ photons released by the jth voxel is captured by the crystal connected by the ith LOR. N represents the total number of voxels to be reconstructed, a il represents the probability that a pair of γ photons released by the lth voxel is captured by the crystal connected by the ith LOR.

[0038] In order to verify the effectiveness of the present application, the Derenzo model is simulated according to the above steps, and the Derenzo model is as shown in Figure 4 . First, the region of interest of the Derenzo model is determined according to step 1, and then the system matrix is calculated by using the conventional solid angle model and the solid angle-distance model of the present application respectively, and the image reconstruction results obtained are as shown in Figure 5 (a) and (b), and the reconstructed image after adding the Gaussian filter is as shown in Figure 5 (c). The calculation of the region of interest and the image reconstruction are performed by using the CUDA architecture for parallel calculation on the GPU, and the calculation time is as shown in Table 1. From Figure 5 and Table 1, it can be seen that the method proposed in the present application not only shortens the image reconstruction time, but also improves the image reconstruction quality, which is beneficial to the accurate judgment of the state of the industrial detection piece.

[0039] Table 1 Time comparison of different methods on GPU

[0040] System matrix Time Traditional solid angle model 350h Region of interest-solid angle model 172h Region of interest-solid angle-distance model 150h

[0041] The above embodiments only exemplarily illustrate the principles and effects of the present application, and part of the applied embodiments, and are not used to limit the present application; it should be noted that, for those skilled in the art, without departing from the concept of the present application, a number of modifications and improvements can be made, which all belong to the protection scope of the present application.

Claims

1. A rapid 3D detection and imaging method for regions of interest in industrial gamma photon detectors, specifically including the following steps: Step 1: Based on the relative position of the object being measured and the detection system, obtain the outline position of the object being measured, i.e., the region of interest, and record and store the voxel information within the region of interest; The contour position of the object under test is obtained based on the position of the workpiece relative to the ring detector. The center of the workpiece is made to coincide with the center of the ring detector. The voxel position information in the contour activity area is calculated. Specifically, the rectangle inscribed in the ring detector is A×B×C, the distance from the center of the workpiece to the contour edge is dx\dy\dz, the size of the reconstructed image is a×b×c, and the distance from the center of the reconstructed image to the contour edge of the object to be reconstructed is tx / ty / tz. The contour position of the workpiece to be reconstructed along the coordinate axes x, y, z is calculated according to equations (1), (2), and (3) to obtain the voxel position information to be reconstructed in the contour activity area. Step 2: Divide the computation system matrix element values ​​into small units that can be executed independently. The calculation of the probability weight corresponding to the voxel in the region of interest for each response line is independent. During image iterative reconstruction, the update of all voxels corresponding to each response line is also independent. Task allocation is performed through the CUDA architecture, assigning a block ID to each response line, and calculating the probability weight corresponding to the voxel in the region of interest in parallel until the computation of the region of interest system matrix is ​​completed. The probability weights corresponding to voxels within the region of interest are calculated as follows: the angle between the voxel center and the detector crystal is the solid angle; the pair of detector crystals corresponding to the response lines passing through the voxel are A and B, respectively; and the minimum solid angle subtended by the detector crystals to the voxel is used as the weight of the system matrix in image reconstruction, i.e., the element value a of the system matrix. ij Then the probability weight of the voxel is μ×a ij In the formula, the weighting coefficient μ is the ratio of the distance between the voxel and detector crystals A and B, μ=D / P; D is the distance between detector crystal A and the center of the voxel corresponding to the smallest solid angle, and P is the distance between the center of the voxel and another detector crystal B; Step 3: Reassign a new block ID to each response line, and perform iterative reconstruction of the 3D-OSEM image based on the region of interest system matrix. The iterative reconstruction formula is as follows: in, S n For the nth subset, y i Let λ represent the i-th LOR. j Let j represent the j-th voxel. and Let a represent the results of the reconstructed image at iterations k+1 and k, respectively. ij Let be an element in the system matrix A, representing the probability that a pair of γ photons emitted by the j-th voxel is captured by the crystal connected to the i-th LOR, and N represent the total number of voxels to be reconstructed. il This represents the probability that a pair of γ photons emitted by the l-th voxel is captured by the crystal connected to the i-th LOR.

2. The method for rapid 3D detection and imaging of regions of interest for industrial gamma photon detectors as described in claim 1, characterized in that: In step 3, during the iterative reconstruction of the image, a Gaussian filter is added during the iterative update.

3. The method for rapid 3D detection and imaging of regions of interest for industrial gamma photon detectors as described in any one of claims 1-2, characterized in that: The CUDA architecture records the resolution of the reconstructed image as M×M, and opens M×M threads in the block. Each block opens a maximum of 1024 threads and a minimum of 32 threads, and the number of threads must be a multiple of 32.

Images