A tomographic gamma scan transmission image reconstruction method combining ART and BM3D

Abstract

This invention discloses a tomographic gamma-ray transmission image reconstruction method combining ART and BM3D, comprising: measuring the transmission gamma energy spectrum and calculating projection data; dividing the tomographic voxels and calculating the voxel attenuation track length; establishing an underdetermined set of equations for transmission image reconstruction; setting initial values ​​for the transmission image; iterating using ART and outputting a noisy image; denoising the noisy image output by ART using BM3D; repeating the iteration and outputting the final transmission image. This solution solves the technical problem of low accuracy in tomographic gamma-ray transmission image reconstruction under conditions of severe undersampling of projection data, achieving transmission image reconstruction while improving its accuracy. This method directly calls the BM3D module in the computer program, making it easy to implement and highly applicable to image reconstruction with undersampling data.

Description

Technical Field

[0001] This invention relates to the field of image data processing and analysis technology, and in particular to a tomographic gamma-ray transmission image reconstruction method combining ART and BM3D. Background Technology

[0002] In nuclear safety monitoring and nuclear waste treatment and disposal, the media, gamma nuclides, and radioactivity levels within waste containers are crucial for the effective identification and classification of radioactive waste containers. According to national standards, gamma radiation detection in radioactive waste containers is a necessary step. Tomographic gamma scanning technology is one of the most advanced non-destructive testing technologies for radioactive waste containers. It uses a gamma detector to scan the waste container at multiple angles, horizontal positions, and heights, reconstructing images of the media attenuation coefficient distribution and nuclide activity distribution within the container, enabling qualitative, quantitative, and locational analysis of gamma nuclides. Tomographic gamma scanning systems are used to inspect 200L nuclear waste containers generated at nuclear fuel manufacturing plants, nuclear waste treatment plants, and nuclear power plants, and are of great significance for the safe handling and disposal of nuclear waste containers.

[0003] Non-destructive testing of nuclear waste containers using tomographic gamma scanning technology is divided into two parts: transmission measurement and emission measurement. In transmission measurement, an external transmission source is used to transmit gamma rays through the waste container, reconstructing a transmission image (i.e., an image of the attenuation coefficient distribution of the medium inside the container), thereby obtaining the attenuation information of gamma rays of different energies within the container. In emission measurement, gamma rays emitted by gamma nuclides inside the container are directly measured, and an emission image (i.e., an image of the nuclide activity distribution within the container) is reconstructed, thereby obtaining the nuclide activity within the container. In emission measurement, gamma rays emitted by nuclides within the container are attenuated and absorbed by the medium inside the container over a distance along the path from the emission location to the detector. Therefore, it is necessary to use the attenuation information of the medium inside the container for gamma rays of different energies from the transmission image to correct the attenuation in the emission image. Thus, the quality of the transmission image has a significant impact on the accuracy of emission image reconstruction and nuclide activity calculation within the container.

[0004] Domestically and internationally, various algorithms are mainly used for tomographic gamma image reconstruction, including Filtered Back Projection (FBP), Expectation-Maximization (EM), Maximum Likelihood Expectation-Maximization (MLEM), Algebraic Reconstruction Algorithm (ART), and Monte Carlo Statistical Iterative Algorithm. Among these, the ART algorithm is the mainstream algorithm for transmission reconstruction. Image spatial resolution and pixel accuracy are two major factors in evaluating the quality of the reconstructed image. Current transmission image reconstruction methods often use large voxels for grid division, resulting in low spatial resolution and poor image quality. To improve image spatial resolution, the image pixel grid needs to be subdivided. However, this results in fewer projected values ​​obtained from transmission measurements than the number of image grids, i.e., undersampling of the projected data. The essence of the transmission image reconstruction problem is solving several linear equations. Under the condition of undersampling of projected data, the transmission image reconstruction equations are severely underdetermined, meaning the number of equations is significantly less than the number of pixels to be solved. Using the traditional ART algorithm to reconstruct transmission images results in low pixel accuracy, rendering traditional iterative algorithms unsuitable.

[0005] To address this issue, a new transmission image reconstruction algorithm is needed that can accurately reconstruct transmission images of finely meshed voxels even under conditions of severe undersampling of projection data in nuclear waste bin tomography gamma-ray transmission measurements. Summary of the Invention

[0006] To address the technical problem of low accuracy in tomographic gamma-scan transmission image reconstruction under severe undersampling of projection data, this invention provides a tomographic gamma-scan transmission image reconstruction method combining ART and BM3D. This method reduces the mean square error of transmission image reconstruction and improves the signal-to-noise ratio of transmission image reconstruction. The method has superior reconstruction effect on transmission images and can still accurately reconstruct transmission images with finely divided voxels even under severe undersampling of transmission measurement projection data. At the same time, this method is easy to implement in a computer.

[0007] This invention is achieved using the following technical solution: A method for reconstructing tomographic gamma-ray transmission images by combining ART and BM3D includes the following steps: Step S1: Based on the transmission source collimator and gamma radiation detector, measure the transmitted gamma energy spectrum and calculate the projection data; Step S2: Divide the fault into voxels and calculate the voxel attenuation track length; Step S3: Define the relationship between projection data, voxel attenuation track length and voxel value in the i-th transmission measurement, and establish an underdetermined set of equations for transmission image reconstruction. Step S4: Set the initial values ​​for the transmission image; Step S5: Iterate based on the ART algorithm and output the noisy image; Step S6: Use BM3D to denoise the noisy image output by the ART algorithm iteration; Step S7: Use the denoised image obtained in step S6 as the initial image value and input it into the ART algorithm iteration process in step S5 again. Repeat steps S5 to S6 to iterate until the preset number of iterations is met, and output the final reconstructed transmission image.

[0008] Specifically, step S1 includes: Step S11: Collimate the γ rays emitted by the γ source using a transmission source collimator to form a narrow beam of γ rays; Step S12: Directly transmit the collimated γ-ray beam through a 200L nuclear waste container, and use a γ-ray radiation detector to measure the γ-ray beam after it penetrates the waste container. One transmission γ-ray energy spectrum is obtained for each transmission measurement. Step S13: Change the horizontal detection position of the gamma radiation detector and the rotation angle of the waste bin, and use the gamma radiation detector... Sub-horizontal movement and waste bins Each rotation, a total of measurements were obtained for one fault. A transmission γ-ray energy spectrum ; Step S14: Remove the waste bin so that the γ-ray beam does not pass through the waste bin and is directly measured by the γ-radiation detector, and obtain an unattenuated γ-ray spectrum that is not blocked by the waste bin. Step S15: Calculate the projection data. The calculation formula is as follows: ; in, E For γ energy, This is the sequence number for each transmission measurement. For the first i Projection data obtained from secondary transmission measurements The area of ​​the characteristic peak in the unattenuated γ-ray spectrum. For the first i The characteristic peak area in the transmission γ-ray energy spectrum obtained from the transmission measurement.

[0009] Specifically, step S2 includes: Step S21: Establish polar coordinates with the fault center of the disc-shaped waste bin fault as the origin; specifically: using The value per unit length is divided along a fault radius of 28cm. A length, =28 / , The unit is cm; The unit angle value is divided into 360° intervals along the fault. From one angle, =360 / The entire fault is divided into A tiny voxel, That is, the total Circles, each circle Individual material blocks; Step S22: Number the voxels sequentially from the outermost circle to the innermost circle in a counter-clockwise direction. The outermost circle is circle 1, denoted as: number 1~ Number; the second circle is represented as: ( +1) number ~ (2× () number; ...; the innermost circle is the number ) Circle, represented as: [( -1)× +1] number~ Number; Step S23: Take the center of the front face of the transmission source collimator as one point and the front face of the gamma radiation detector as another point, connect the two points to form a straight line, and calculate the length of this straight line passing through each voxel, which is the attenuation track length of the transmitted gamma ray beam in each voxel.

[0010] Specifically, the relationship in step S3 is represented as follows: ; in, For the first Projection data obtained from secondary transmission measurements For the first In the second transmission measurement, the X-ray beam passes through the first... Attenuation track length at voxel number 1 , Transmission γ-ray energy spectrum quantity Much smaller than the number of tomographic voxels , For the first The pixel value of a voxel; The underdetermined equations for the reconstruction of the transmission image are constructed in conjunction with the total number of lines passed through a tomography section. Secondary transmission measurement, and obtain The equation for the projected data is as follows: ; ; ; … .

[0011] Specifically, the initial values ​​of the transmission image in step S4 are set as follows: [ , , …, = [0.001, 0.001, …, 0.001]; in, , , …, They are respectively No. 1, No. 2, ..., No. The initial pixel value of the tomographic voxel.

[0012] Specifically, the iteration based on the ART algorithm in step S5 includes: Step S51: Input the initial values ​​of the transmission image into the ART algorithm, for the first... For the transmission measurement, the pixel value of the j-th voxel after the k-th iteration is calculated using the following formula: ; in, Let be the pixel value of the j-th voxel after the k-th iteration. Let $j$ be the pixel value of the $j$-th voxel after the $k-1$-th iteration. When $k=1$, Let be the initial pixel value of the j-th voxel, and λ be a relaxation factor controlling the iteration step size. The projection data obtained from the i-th transmission measurement. For the first The attenuation track length of the X-ray beam when it passes through voxel j in the secondary transmission measurement; For each transmission measurement, the ART algorithm is used to iterate once using the above formula, for a total of Sub-transmission measurement, when the iteration number k= One iteration is completed in time; Step S52: After passing through the ART algorithm After round of iterations, the output is... The pixel values ​​of the voxels are as follows: , , …, , which is the noisy image Z.

[0013] Specifically, step S6 includes: Step S61: Extract the noisy image from the ART algorithm output. The input image is used for block matching; where the reference block has a side length of... The search window has a side length of 1. The sliding step size is ; Step S62: Process the three-dimensional matrix obtained in step S61 Perform collaborative filtering; Step S63: Estimate the denoised 3D matrix obtained in step S62. As input, an aggregation operation is performed to return all the denoised image patches to the original noisy image. Corresponding spatial location The weighted average of the estimates for overlapping blocks is calculated using the following formula: ; in, For the fault number The first circle The position coordinates of the voxels , ; Based on the denoised image in the corresponding The estimated pixel value of the location. The estimated value of the denoised image; The denoised 3D matrix is ​​then restored to its corresponding... Pixel estimate of the location; To be assigned Aggregate weights; Step S64: The estimated value of the base denoised image obtained in step S63. Using the guiding input image, the block matching operation in step S61 is repeated to obtain the basic denoising 3D matrix. And based on the aforementioned basic denoising three-dimensional matrix The coordinates of stacked image patches in the original noisy image Image patches are extracted and stacked to obtain a noisy 3D matrix. ; Step S65: Apply the basic denoised 3D matrix obtained in step S64 and noisy three-dimensional matrix Perform Wiener filtering operation; Step S66: The final denoised 3D matrix estimate obtained in step S65 is... As input, the aggregation operation in step S63 is repeated, and the final image after denoising by the BM3D algorithm is output. .

[0014] Specifically, the block matching operation in step S61 further includes the following sub-steps: Step S611: Process the ART algorithm... The size output after round iteration is Noisy images Divided into a series of sizes Reference block ; Step S612: Set the size to Within the search window, by step size Slide selection and Similar image patches And by normalizing the Euclidean distance The formula for measuring their similarity is: ; in, The total number of pixels in the two-dimensional reference block and , The smaller the value, the more similar the two image patches are; Step S613: Set the similarity threshold Retain those with similarity greater than this threshold, i.e. Blocks smaller than this threshold are selected, and all similar blocks are stacked along the depth direction to form a block of size [size missing]. 3D matrix ,in The number of similar blocks to retain.

[0015] Specifically, step S62 includes the following sub-steps: Step S621: Through three-dimensional orthogonal transformation Will Transform to the frequency domain to obtain the transform coefficients. Hard threshold shrinkage operator is used Directly truncate the transform coefficients The noise in is represented as The specific form is as follows: ; in, These are the transform coefficients after denoising. The standard deviation of noise. This is the hard threshold parameter; Step S622: Analyze the denoised transform coefficients. Perform three-dimensional inverse transformation Obtain the denoised 3D matrix estimate. .

[0016] Specifically, step S65, the Wiener filtering operation, includes the following sub-steps: Step S651: Through three-dimensional orthogonal transformation Will and Transform to the frequency domain to obtain the corresponding fundamental transform coefficients. and noisy transform coefficients Then Wiener filtering is performed, with the following formula: ; in, These are the final transform coefficients after Wiener filtering; The noise variance is estimated based on the noise standard deviation σ. Step S652: For the final transformation coefficients Perform three-dimensional inverse transformation The final denoised 3D matrix estimate is obtained. .

[0017] The beneficial effects of this invention are as follows: This solution solves the technical problem of low accuracy in tomographic gamma-scan transmission image reconstruction under severe undersampling of projection data by designing a tomographic gamma-scan transmission image reconstruction method that combines ART and BM3D. It accurately and effectively realizes transmission image reconstruction and improves the accuracy of transmission image reconstruction. At the same time, this method directly calls the BM3D module in the computer program, which is easy to implement and has promotional value and strong universality for image reconstruction with undersampling data. Attached Figure Description

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

[0019] Figure 1 This is a basic flowchart of the tomographic gamma-scan transmission image reconstruction method in an embodiment of the present invention; Figure 2 This is a schematic diagram of a tomographic sample in an embodiment of the present invention; Figure 3 This is a reference image (0.662 MeV) of the tomographic sample in an embodiment of the present invention. Figure 4 This is a reference image (1.173 MeV) of the tomographic sample in an embodiment of the present invention. Figure 5 This is a reference image (1.332 MeV) of the tomographic sample in an embodiment of the present invention. Figure 6 This is a transmission image (0.662 MeV) reconstructed by ART in an embodiment of the present invention. Figure 7 This is a transmission image (1.173 MeV) reconstructed by ART in an embodiment of the present invention. Figure 8 This is a transmission image (1.332 MeV) reconstructed by ART in an embodiment of the present invention. Figure 9 This is a transmission image (0.662 MeV) jointly reconstructed by ART and BM3D in an embodiment of the present invention. Figure 10This is a transmission image (1.173 MeV) jointly reconstructed by ART and BM3D in an embodiment of the present invention. Figure 11 This is a transmission image (1.332 MeV) jointly reconstructed by ART and BM3D in an embodiment of the present invention.

[0020] Wherein: 1-Outer edge of waste bin sample, A-Polyethylene, B-Polyvinyl chloride, C-Concrete, D-Air. Detailed Implementation

[0021] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. The components of the embodiments of the present invention described and shown in the accompanying drawings can generally be arranged and designed in various different configurations.

[0022] It should be noted that similar labels and letters in the following figures indicate similar items. Therefore, once an item is defined in one figure, it does not need to be further defined and explained in subsequent figures.

[0023] The following is in conjunction with the appendix Figure 1 ~Attached Figure 11 The following describes some embodiments of the present invention in detail. Unless otherwise specified, the following embodiments and features can be combined with each other.

[0024] This invention proposes a method for reconstructing tomographic gamma-ray transmission images by combining ART and BM3D. In a preferred embodiment, the method includes the following steps: Step S1: Based on the transmission source collimator and gamma radiation detector, measure the transmitted gamma energy spectrum and calculate the projection data; Step S2: Divide the fault into voxels and calculate the voxel attenuation track length; Step S3: Define the relationship between projection data, voxel attenuation track length and voxel value in the i-th transmission measurement, and establish an underdetermined set of equations for transmission image reconstruction. Step S4: Set the initial values ​​for the transmission image; Step S5: Iterate based on the ART algorithm and output the noisy image; Step S6: Use BM3D to denoise the noisy image output by the ART algorithm iteration; Step S7: Use the denoised image obtained in step S6 as the initial image value and input it into the ART algorithm iteration process in step S5 again. Repeat steps S5 to S6 to iterate until the preset number of iterations is met, and output the final reconstructed transmission image.

[0025] like Figure 1 As shown, the basic steps of this method include: S1: Measure the transmitted γ-ray energy spectrum and calculate the projection data; S2: Divide the fault voxels and calculate the voxel attenuation track length; S3: Establish the underdetermined system of equations for reconstructing the transmission image; S4: Set the initial values ​​for the transmission image; S5: Iterate using ART and output a noisy image; S6: Use BM3D to denoise the noisy image output by ART; S7: Repeat the iteration and output the final transmission image.

[0026] Using the above methods, when the fault is divided into small grid voxels, the number of measured transmission γ-ray spectra is much smaller than the number of fault voxels, which is a sparse sampling method that can quickly complete the transmission measurement of the waste bin; under the condition of severe undersampling of the transmission measurement projection data, accurate reconstruction of the transmission image can be achieved.

[0027] The steps are explained in detail below with reference to specific embodiments: S1: Measure the transmitted gamma spectrum and calculate the projection data, specifically including: S11: Using a transmission source collimator to target a γ source with high radioactivity. and The emitted gamma rays are collimated to form a narrow gamma-ray beam. The activity of the source is 3.69 × 10⁻⁶. 7 Bq emits gamma rays with an energy of 0.662 MeV. The activity of the source is 3.64 × 10⁻⁶. 7 Bq emits gamma rays with energies of 1.173 MeV and 1.332 MeV; S12: A tomographic sample of a nuclear waste bin was constructed using three materials of different densities: polyethylene, polyvinyl chloride, and concrete. Then, a collimated gamma-ray beam was directly transmitted through the sample. The interior of the sample was as follows: Figure 2 As shown, the diameter of polyethylene is 12cm, the diameter of polyvinyl chloride is 10cm, the diameter of concrete is 9cm, and the diameter of the fracture sample is 56cm. The gamma-ray beam after penetrating the waste can sample is measured using a 2-inch sodium iodide detector, and one transmission gamma spectrum is obtained from one transmission measurement. S13: By changing the horizontal detection position of the sodium iodide detector and the rotation angle of the waste bin, and through 4 horizontal movements of the sodium iodide detector and 24 rotations of the waste bin, a total of 96 transmission gamma spectra were obtained by one tomographic measurement, 96=4×24. S14: Remove the waste bin so that the gamma-ray beam does not pass through the waste bin and is directly measured by the gamma radiation detector, and obtain an unattenuated gamma energy spectrum that is not blocked by the waste bin. S15: Calculate the projection data, using the following formula: ; in, For γ energy, This is the sequence number for each transmission measurement. For the first Projection data obtained from secondary transmission measurements The area of ​​the characteristic peak in the unattenuated γ-ray spectrum. For the first The characteristic peak area in the transmission γ energy spectrum obtained from the transmission measurement, combined with the 96 transmission γ energy spectra obtained in step S13, yields 96 projection data for each γ energy at the three energies of 0.662 MeV, 1.173 MeV, and 1.332 MeV.

[0028] S2: Divide the fault into voxels and calculate the voxel attenuation track length, specifically including: S21: The fracture of the waste bin sample is disc-shaped. Polar coordinates are established with the center of the fracture as the origin. The length value is 2.333cm. The 28cm radius of the fracture can be divided into 12 lengths. The angle value is 15°. The 360° of the fracture can be divided into 72 angles. The entire fracture is divided into 864 small voxels. 864 = 12 × 72, that is, a total of 12 rings, with 72 voxel blocks in each ring. S22: Number the voxels in a counter-clockwise direction from the outer circle to the inner circle. That is, the first circle (outermost circle) is numbered 1 to 72, the second circle is numbered 73 to 144, ..., and the 12th circle (innermost circle) is numbered 793 to 864. S23: Take the center of the front face of the transmission source collimator as one point and the front face of the sodium iodide detector as another point. Connect these two points to form a straight line. Calculate the length of this straight line passing through each voxel, which is the attenuation track length of the transmitted γ-ray beam in each voxel.

[0029] S3: Establish the underdetermined equations for reconstructing the transmission image, specifically including: No. In the secondary transmission measurement, the relationship between the projection data, voxel attenuation track length, and voxel value is as follows: ; in, For the first Projection data obtained from secondary transmission measurements For the first In the second transmission measurement, the X-ray beam passes through the first... The attenuation track length at voxel number 1 ≤ ≤96, 1≤ The number of transmitted gamma spectra measured (≤864) is much smaller than the number of tomographic voxels (864), therefore it is sparse sampling. For the first The pixel value of a voxel.

[0030] Combining 96 transmission measurements and 96 projection data points obtained from a single fault, the underdetermined equations for reconstructing a transmission image with an energy of 0.662 MeV are as follows: ; ; ; … ; The underdetermined equations for reconstructing a transmission image with an energy of 1.173 MeV are as follows: ; ; ; … ; The underdetermined equations for reconstructing a transmission image with an energy of 1.332 MeV are as follows: ; ; ; … .

[0031] S4: Set the initial values ​​for the transmission image, specifically including: When E is 0.662 MeV, 1.173 MeV, or 1.332 MeV, the initial values ​​for the transmitted image are set as follows: [ , , …, = [0.001, 0.001, …, 0.001]; in, , , …, These are the initial pixel values ​​for voxels numbered 1, 2, ..., 864.

[0032] S5: Iterate using the ART algorithm and output a noisy image, specifically including: S51: Input the initial values ​​of the transmission image in step S4 into the ART algorithm. For each transmission measurement, use the ART algorithm to perform one iteration. One round of iteration is completed through 96 calculations. The relaxation factor λ is 0.05. S52: After 5 iterations of the ART algorithm, the pixel values ​​of the 864 voxels are as follows: , ,…, This is a noisy image. .

[0033] S6: Denoising the noisy image output by ART using BM3D, specifically including: S61: The noisy image output by the ART algorithm As input, the image is used for block matching, and a 3D matrix is ​​output. ; where the reference block side length =8, meaning the reference block size is 8×8; search window side length =11, meaning the search window size is 11×11; sliding step size =3; Number of similar blocks to retain =16; S62: The three-dimensional matrix obtained in step S61 Perform collaborative filtering to obtain the denoised 3D matrix estimate. Among them, the hard threshold parameter =2.5; Noise standard deviation The initial value is 0.05; S63: The denoised 3D matrix estimate obtained in step S62 As input, perform aggregation operations and output an estimate of the base denoised image. ; S64: Estimate the base denoised image obtained in S63 Using the guiding input image, the block matching operation in step S61 is repeated to obtain the basic denoising 3D matrix. and according to The coordinates of stacked image patches in the original noisy image Image patches are extracted and stacked to obtain a noisy 3D matrix. ; S65: Apply the basic denoised 3D matrix obtained in step S64. and noisy three-dimensional matrix Perform Wiener filtering to obtain the final denoised 3D matrix estimate. ; S66: The final denoised 3D matrix estimate obtained in step S65 is... As input, the aggregation operation in step S63 is repeated, and the final image after denoising by the BM3D algorithm is output. .

[0034] S7: Repeat the iteration and output the final transmission image, specifically including: The denoised image obtained in step S6 is used as the initial image value and fed into the ART iteration process for the next iteration. Steps S5 and S6 are repeated until the total number of iterations reaches the preset number (e.g., 100 times), and the final reconstructed transmission images with energies of 0.662 MeV, 1.173 MeV, and 1.332 MeV are output.

[0035] In this embodiment, a cross-sectional sample of a waste bin is constructed, utilizing... and Reference images of energies of 0.662 MeV, 1.173 MeV, and 1.332 MeV were obtained from transmission measurements of the source, as shown in the image. Figures 3-5 As shown.

[0036] Using the ART algorithm, transmission images with energies of 0.662 MeV, 1.173 MeV, and 1.332 MeV were reconstructed, such as... Figures 6-8 As can be seen, the single ART algorithm cannot accurately reconstruct the transmission image, and there are many artifacts in the image, resulting in an unclear image.

[0037] Using a combined ART and BM3D method, transmission images with energies of 0.662 MeV, 1.173 MeV, and 1.332 MeV were reconstructed, such as... Figures 9-11 As can be seen, the method of combining ART and BM3D can accurately reconstruct transmission images with fewer artifacts and higher transmission image quality.

[0038] Meanwhile, the mean square error (MSE) and signal-to-noise ratio (SNR) are used as evaluation indicators for the quality of transmission image reconstruction, as shown in the following formula: ; ; in, The first in the transmission image The pixel reconstruction value of the individual pixels, For the reference image, the first The pixel value of a voxel.

[0039] The MSE of the reconstructed 0.662MeV, 1.173MeV, and 1.332MeV energy transmission images using the joint ART and BM3D method were 0.64, 0.60, and 0.64 times that of the ART algorithm, respectively, reducing the mean square error of the transmission image reconstruction. The SNR of the reconstructed 0.662MeV, 1.173MeV, and 1.332MeV energy transmission images using the joint ART and BM3D method were 1.50, 1.56, and 1.48 times that of the ART algorithm, respectively, improving the signal-to-noise ratio of the transmission image reconstruction.

[0040] In summary, both qualitative and quantitative analyses of transmission image reconstruction demonstrate that the quality of transmission image reconstruction using the combined ART and BM3D method is higher than that of the ART algorithm alone.

[0041] For the foregoing embodiments, in order to simplify the description, they are all described as a series of actions. However, those skilled in the art should understand that this application is not limited to the described order of actions, because according to this application, some steps can be performed in other orders or simultaneously. Furthermore, those skilled in the art should also understand that the embodiments described in the specification are preferred embodiments, and the actions involved are not necessarily essential to this application.

[0042] The above embodiments describe the basic principles, main features, and advantages of the present invention. Those skilled in the art should understand that the present invention is not limited to the above embodiments. The embodiments and descriptions in the specification are merely illustrative of the principles of the invention. Modifications and variations made by those skilled in the art without departing from the spirit and scope of the invention should be within the protection scope of the appended claims.

Claims

1. A method for reconstructing tomographic gamma-ray transmission images combining ART and BM3D, characterized in that, Includes the following steps: Step S1: Based on the transmission source collimator and gamma radiation detector, measure the transmitted gamma energy spectrum and calculate the projection data; Step S2: Divide the fault into voxels and calculate the voxel attenuation track length; Step S3: Define the relationship between projection data, voxel attenuation track length and voxel value in the i-th transmission measurement, and establish an underdetermined set of equations for transmission image reconstruction. Step S4: Set the initial values ​​for the transmission image; Step S5: Iterate based on the ART algorithm and output the noisy image; Step S6: Use BM3D to denoise the noisy image output by the ART algorithm iteration; Step S7: Use the denoised image obtained in step S6 as the initial image value and input it into the ART algorithm iteration process in step S5 again. Repeat steps S5 to S6 to iterate until the preset number of iterations is met, and output the final reconstructed transmission image.

2. The method for reconstructing tomographic gamma-ray transmission images combining ART and BM3D as described in claim 1, characterized in that, Step S1 specifically includes: Step S11: Collimate the γ rays emitted by the γ source using a transmission source collimator to form a narrow beam of γ rays; Step S12: Directly transmit the collimated γ-ray beam through a 200L nuclear waste container, and use a γ-ray radiation detector to measure the γ-ray beam after it penetrates the waste container. One transmission γ-ray energy spectrum is obtained for each transmission measurement. Step S13: Change the horizontal detection position of the gamma radiation detector and the rotation angle of the waste bin, and use the gamma radiation detector... Sub-horizontal movement and waste bins Each rotation, a total of measurements were obtained for one fault. A transmission γ-ray energy spectrum ; Step S14: Remove the waste bin so that the γ-ray beam does not pass through the waste bin and is directly measured by the γ-radiation detector, and obtain an unattenuated γ-ray spectrum that is not blocked by the waste bin. Step S15: Calculate the projection data. The calculation formula is as follows: ； in, E For γ energy, This is the sequence number for each transmission measurement. For the first i Projection data obtained from secondary transmission measurements The area of ​​the characteristic peak in the unattenuated γ-ray spectrum. For the first i The characteristic peak area in the transmission γ-ray energy spectrum obtained from the transmission measurement.

3. The method for reconstructing tomographic gamma-ray transmission images combining ART and BM3D as described in claim 2, characterized in that, Step S2 specifically includes: Step S21: Establish polar coordinates with the fault center of the disc-shaped waste bin fault as the origin; specifically: using The value per unit length is divided along a fault radius of 28cm. A length, =28 / , The unit is cm; The unit angle value is divided into 360° intervals along the fault. From one angle, =360 / The entire fault is divided into A tiny voxel, That is, the total Circles, each circle Individual material blocks; Step S22: Number the voxels sequentially from the outermost circle to the innermost circle in a counter-clockwise direction. The outermost circle is circle 1, denoted as: number 1~ Number; the second circle is represented as: ( +1) number ~ (2× () number; ...; the innermost circle is the number ) Circle, represented as: [( -1)× +1] number~ Number; Step S23: Take the center of the front face of the transmission source collimator as one point and the front face of the gamma radiation detector as another point, connect the two points to form a straight line, and calculate the length of this straight line passing through each voxel, which is the attenuation track length of the transmitted gamma ray beam in each voxel.

4. The method for reconstructing tomographic gamma-ray transmission images combining ART and BM3D as described in claim 3, characterized in that, The relationship in step S3 is represented as follows: ； in, For the first Projection data obtained from secondary transmission measurements For the first In the second transmission measurement, the X-ray beam passes through the first... Attenuation track length at voxel number 1 , Number of transmission γ-ray spectra Much smaller than the number of tomographic voxels , For the first The pixel value of a voxel; The underdetermined equations for the reconstruction of the transmission image are constructed in conjunction with the total number of lines passed through a tomography section. Secondary transmission measurement, and obtain The equation for the projected data is as follows: ； ； ； … 。 5. The method for reconstructing tomographic gamma-ray transmission images combining ART and BM3D as described in claim 4, characterized in that, The initial values ​​of the transmission image in step S4 are set as follows: [ , , …, ]=[0.001, 0.001, …,0.001]； in, , , …, They are respectively No. 1, No. 2, ..., No. The initial pixel value of the tomographic voxel.

6. The method for reconstructing tomographic gamma-ray transmission images combining ART and BM3D as described in claim 5, characterized in that, The iteration based on the ART algorithm in step S5 specifically includes: Step S51: Input the initial values ​​of the transmission image into the ART algorithm, for the first... For the transmission measurement, the pixel value of the j-th voxel after the k-th iteration is calculated using the following formula: ； in, Let be the pixel value of the j-th voxel after the k-th iteration. Let $j$ be the pixel value of the $j$-th voxel after the $k-1$-th iteration. When $k=1$, Let be the initial pixel value of the j-th voxel, and λ be a relaxation factor controlling the iteration step size. The projection data obtained from the i-th transmission measurement. For the first The attenuation track length of the X-ray beam when it passes through voxel j in the secondary transmission measurement; For each transmission measurement, the ART algorithm is used to iterate once using the above formula, for a total of Sub-transmission measurement, when the iteration number k= One iteration is completed in time; Step S52: After passing through the ART algorithm After round of iterations, the output is... The pixel values ​​of the voxels are as follows: , , …, , which is the noisy image Z.

7. The method for reconstructing tomographic gamma-ray transmission images combining ART and BM3D as described in claim 6, characterized in that, Step S6 specifically includes: Step S61: Extract the noisy image from the ART algorithm output. The input image is used for block matching; where the reference block has a side length of... The search window has a side length of 1. The sliding step size is ; Step S62: Process the three-dimensional matrix obtained in step S61 Perform collaborative filtering operation; Step S63: Estimate the denoised 3D matrix obtained in step S62. As input, an aggregation operation is performed to return all the denoised image patches to the original noisy image. Corresponding spatial location The weighted average of the estimates for overlapping blocks is calculated using the following formula: ； in, For the fault number The first circle The position coordinates of the voxels , ; Based on the denoised image in the corresponding The estimated pixel value of the location. The estimated value of the denoised image as a whole; The denoised 3D matrix is ​​then restored to its corresponding... Pixel estimate of the location; To be assigned Aggregate weights; Step S64: The estimated value of the base denoised image obtained in step S63. Using the guiding input image, the block matching operation in step S61 is repeated to obtain the basic denoising 3D matrix. And based on the aforementioned basic denoising three-dimensional matrix The coordinates of stacked image patches in the original noisy image Image patches are extracted and stacked to obtain a noisy 3D matrix. ; Step S65: Apply the basic denoised 3D matrix obtained in step S64 and noisy three-dimensional matrix Perform Wiener filtering operation; Step S66: The final denoised 3D matrix estimate obtained in step S65 is... As input, the aggregation operation in step S63 is repeated, and the final image after denoising by the BM3D algorithm is output. .

8. The method for reconstructing tomographic gamma-ray transmission images combining ART and BM3D as described in claim 7, characterized in that, The block matching operation in step S61 further includes the following sub-steps: Step S611: Process the ART algorithm... The size output after round iteration is Noisy images Divided into a series of sizes Reference block ; Step S612: Set the size to Within the search window, by step size Slide selection and Similar image patches And by normalizing the Euclidean distance The formula for measuring their similarity is: ； in, The total number of pixels in the two-dimensional reference block and , The smaller the value, the more similar the two image patches are; Step S613: Set the similarity threshold Retain those with similarity greater than this threshold, i.e. Blocks smaller than this threshold are selected, and all similar blocks are stacked along the depth direction to form a block of size [size missing]. 3D matrix ,in The number of similar blocks to retain.

9. The method for reconstructing tomographic gamma-ray transmission images combining ART and BM3D as described in claim 8, characterized in that, Step S62 includes the following sub-steps: Step S621: Through three-dimensional orthogonal transformation Will Transform to the frequency domain to obtain the transform coefficients. Hard threshold shrinkage operator is used Directly truncate the transform coefficients The noise in is represented as The specific form is as follows: ； in, These are the transform coefficients after denoising. The standard deviation of noise. This is the hard threshold parameter; Step S622: Analyze the denoised transform coefficients. Perform three-dimensional inverse transformation Obtain the denoised 3D matrix estimate. .

10. The method for reconstructing tomographic gamma-ray transmission images combining ART and BM3D as described in claim 9, characterized in that, The Wiener filtering operation in step S65 includes the following sub-steps: Step S651: Through three-dimensional orthogonal transformation Will and Transform to the frequency domain to obtain the corresponding fundamental transform coefficients. and noisy transform coefficients Then Wiener filtering is performed, with the following formula: ； in, These are the final transform coefficients after Wiener filtering; The noise variance is estimated based on the noise standard deviation σ. Step S652: For the final transformation coefficients Perform three-dimensional inverse transformation The final denoised 3D matrix estimate is obtained. .

Images