An Industrial CT Three-Dimensional Image Reconstruction Method Based on Data Rearrangement and Conjugate Rays

Abstract

The present invention discloses an industrial CT three-dimensional image reconstruction method based on data rearrangement and conjugate rays. The method of the present invention first obtains the projection data of the sample to be measured through circular orbit scanning of industrial cone-beam CT; then, the projection data is rearranged by converting the planar detector into a cylindrical detector, and four-neighbor bilinear interpolation is used to improve the accuracy of the projection data; then, the rearranged projection data is cosine weighted, and a three-dimensional empirical weighting function is constructed by a method based on conjugate rays for secondary weighting; then, the weighted projection data is filtered by a one-dimensional ramp filter, and the reconstruction result is obtained through back-projection. The present invention has strong robustness, can eliminate the reconstruction artifacts of projection data with a large cone angle, and thus improve the accuracy of industrial CT measurement.

Description

Technical Field

[0001] The present invention relates to an industrial CT three-dimensional image reconstruction method based on data rearrangement and conjugate rays, and more particularly to an approximate reconstruction method for industrial cone-beam CT circular orbit three-dimensional images. Background Art

[0002] Industrial CT is the industrial application of computed tomography technology. High-energy X-rays are used to penetrate the workpiece to be detected to obtain projection data, and then two-dimensional or three-dimensional images of the workpiece are obtained through a reconstruction algorithm. Industrial CT can clearly, accurately and intuitively display the size, position, shape, composition, material and defect conditions of the internal structure of the workpiece to be detected without contact and damage, and is not affected by the material, shape and surface of the workpiece. Industrial CT has the advantages of fast detection speed, high spatial and density resolution, etc. Therefore, industrial CT is widely used in the field of industrial non-destructive testing.

[0003] Currently, industrial CT mainly uses cone-beam X-rays. During scanning, the scanning path of the X-ray source is a circular orbit relative to the workpiece to be detected. This scanning mode has the advantages of high scanning speed and radiation utilization rate, can simplify the hardware structure of industrial CT, and improve the reconstruction speed at the same time. However, circular orbit scanning cannot meet the Tuy condition for cone-beam CT reconstruction, so accurate reconstruction of three-dimensional images cannot be achieved, and only approximate reconstruction can be realized through algorithms.

[0004] For industrial cone-beam CT circular orbit scanning, there are currently improved three-dimensional Radon inverse transform, improved Grangeat algorithm and FDK algorithm. Among them, both the improved three-dimensional Radon inverse transform and the improved Grangeat algorithm were applicable to the accurate reconstruction of three-dimensional images before improvement, and after improvement, the algorithms were made applicable to approximate reconstruction by interpolation. The processing methods of these two algorithms for projection data are relatively cumbersome and complex, and interpolation is required for the missing data of the Radon transform; the FDK algorithm is accurate for the reconstruction of the central plane of the workpiece to be detected, but as the X-ray cone angle increases, the error of the reconstructed result of the penetrated plane will gradually increase, and finally obvious artifacts will occur. Summary of the Invention

[0005] The present invention proposes an industrial CT three-dimensional image reconstruction method based on data rearrangement and conjugate rays in view of the deficiencies of existing industrial CT reconstruction algorithms. This method first obtains the projection data of the measured sample through circular orbit scanning of industrial cone-beam CT; then, the projection data is rearranged by converting the planar detector into a cylindrical detector, and four-neighbor bilinear interpolation is used to improve the accuracy of the projection data; then, the rearranged projection data is cosine-weighted, and a three-dimensional empirical weighting function is constructed by a method based on conjugate rays for secondary weighting; then, a one-dimensional ramp filter is used to filter the weighted projection data, and the reconstruction result is obtained through back-projection. The present invention has strong robustness and can eliminate the reconstruction artifacts of projection data with a large cone angle, thereby improving the accuracy of industrial CT measurement.

[0006] The technical solution adopted by the present invention is an industrial CT three-dimensional image reconstruction method based on data rearrangement and conjugate rays, which is specifically implemented according to the following steps:

[0007] Step 1: Detect the workpiece by using the circular orbit scanning mode for the cone-beam X-rays of the industrial CT to obtain the projection data of the planar detector at different scanning angles;

[0008] Step 2: Perform a homothetic transformation on the planar detector in Step 1 to obtain a virtual planar detector smaller than the real detector, transfer the projection data on the real detector to the virtual planar detector, and then convert the three-dimensional cone-beam X-rays in Step 1 into multiple two-dimensional virtual cone-beam X-ray sources, and further convert the virtual planar detector into a virtual cylindrical detector;

[0009] Step 3: Use four-neighbor bilinear interpolation for the projection data on the virtual planar detector in Step 2, and fill the interpolated projection data into the virtual cylindrical detector in Step 2;

[0010] Step 4: Perform pre-weighting on the projection data obtained in Step 3 by using the cosine function of the X-ray incident angle;

[0011] Step 5: Construct a three-dimensional empirical weighting function by a method based on conjugate rays to perform secondary weighting on the data pre-weighted in Step 4;

[0012] Step 6: Use the fast Fourier transform to implement the convolution operation of the one-dimensional ramp filter and the data after secondary weighting in Step 5, so as to achieve the filtering effect;

[0013] Step 7: Perform a back-projection operation on the filtered data in Step 6 to obtain the value of the reconstruction point;

[0014] Step 8: Perform operations on all points to be reconstructed according to Steps 4-Step 7, and finally obtain the three-dimensional reconstruction result.

[0015] The beneficial effects of the present invention are as follows: By adopting the method of data rearrangement, the operation process of reconstruction is simplified. Using the theory and empirical method of conjugate rays, different contribution values are set for each X-ray passing through the reconstruction point, and finally the reconstruction of the industrial CT three-dimensional image is realized. The method of the present invention has a relatively high reconstruction efficiency, can eliminate the reconstruction artifacts with a large cone angle, and improves the accuracy of industrial CT measurement; at the same time, the method of the present invention has strong robustness and can be adapted to the reconstruction of projections of workpieces with various shapes. Brief Description of the Drawings

[0016] Figure 1 is the flowchart of the steps of the method of the present invention.

[0017] Figure 2 is the schematic diagram of industrial CT scanning of the method of the present invention.

[0018] Figure 3 is the schematic diagram of the virtual planar detector of the method of the present invention.

[0019] Figure 4 is the schematic diagram of the detector plane coordinate system after the rotary table rotates in the method of the present invention.

[0020] Figure 5 is the schematic diagram of the virtual cylindrical detector of the method of the present invention.

[0021] Figure 6 is the schematic diagram of the interpolation of projection data of the virtual screen detector of the method of the present invention.

[0022] Figure 7 is the schematic diagram of the conjugate ray of the method of the present invention.

[0023] Figure 8 is the physical diagram of the standard sphere plate of the method of the present invention.

[0024] Figure 9 is the three-dimensional reconstruction diagram of the standard sphere plate after industrial CT scanning of the method of the present invention. Detailed Embodiment

[0025] The present invention will be further described below with reference to the accompanying drawings.

[0026] As Figure 1 shown, the steps of the method of the present invention are as follows:

[0027] Step 1: Detect the workpiece by using the circular orbit scanning mode for the cone-beam X-rays of the industrial CT to obtain the projection data of the planar detector at different scanning angles;

[0028] As Figure 2As shown, the X-ray source emits cone-beam X-rays that penetrate the workpiece to be detected, and the flat-panel detector senses the attenuation of the X-ray intensity to obtain a projection data map. Then, the rotating table rotates clockwise. For each 1° rotation, the flat-panel detector can acquire 1 projection data map. During the measurement process, when the rotating table rotates one full circle, the flat-panel detector can obtain a total of 360 projection data maps, and the cone-beam X-rays form a circular scanning orbit relative to the rotating table. At the same time, with Figure 1 the coordinate system shown, a scanning object coordinate system is established.

[0029] Step 2: Perform a homothetic transformation on the planar detector in Step 1 to obtain a virtual planar detector smaller than the real detector, transfer the projection data on the real detector to the virtual planar detector, then convert the three-dimensional cone-beam X-rays in Step 1 into multiple two-dimensional virtual cone-beam X-ray sources, and further convert the virtual planar detector into a virtual cylindrical detector;

[0030] As Figure 3 shown, with the X-ray source S as the homothetic center, perform a homothetic transformation on the real planar detector so that it coincides with the central plane of the rotating table in Step 1 to obtain a virtual planar detector. Establish a space rectangular coordinate system with the center O of the virtual planar detector, and transfer the projection data on the real detector to the virtual planar detector. As Figure 4 shown, after the rotating table rotates clockwise, let a be the horizontal direction of the virtual planar detector, b be the direction of the line connecting the X-ray source S and the center O of the virtual planar detector, c be the vertical direction of the virtual planar detector, which coincides with the z direction of the scanning object coordinate system, be the angle of counterclockwise rotation of the virtual planar detector relative to the rotating table, and this angle is equal to the angle of clockwise rotation of the rotating table. Then, the conversion relationship between the virtual planar detector coordinate system and the scanning object is

[0031]

[0032] Then, establish the conversion relationship between the virtual planar detector parameter space and the virtual cylindrical detector parameter space (s, c, θ)

[0033]

[0034] In the formula, R is the distance from the X-ray source S to the center O of the virtual planar detector. As Figure 5As shown, this coordinate transformation rearranges the projection points on the virtual planar detector to form a curve on the same horizontal plane, and then connects the curves to form a surface, finally converting the virtual planar detector into a virtual cylindrical detector. In the virtual cylindrical detector, s is the direction perpendicular to the line connecting the X-ray source S' and the center O of the virtual cylindrical detector, t is the direction of the line connecting the X-ray source S' and the center O of the virtual cylindrical detector, and θ is the angle of counterclockwise rotation of the s direction relative to the x direction of the scanning object coordinate system.

[0035] Step 3: Perform four-neighbor bilinear interpolation on the projection data on the virtual planar detector in Step 2, and then fill the virtual cylindrical detector in Step 2 with the interpolated projection data;

[0036] As Figure 6 shown, P is the projection on the virtual planar detector, and let its projection data be The projection data of the detector unit adjacent to the lower left of this point is The projection data of the detector unit adjacent to the upper left is The projection data of the detector unit adjacent to the lower right is The projection data of the detector unit adjacent to the upper right is Calculate the coordinate difference between p and p 11 of

[0037] u = a - a0 (5)

[0038] v = c - c0 (6)

[0039] Then, perform four-neighbor bilinear interpolation on and there is

[0040]

[0041] Fill the virtual cylindrical detector in Step 2 with the interpolated projection data, and make the number of filling points of the curve at each height be an integer power of 2. The filling relationship is

[0042]

[0043] where is the projection data after filling the virtual cylindrical detector.

[0044] Step 4: Perform pre-weighting on the projection data obtained in Step 3 using the cosine function of the X-ray incident angle;

[0045] Let the coordinates of the reconstruction point of the workpiece to be detected be (x, y, z), and the value of h in the coordinates of the projection point on the virtual cylindrical detector can be calculated from the triangle similarity relationship

[0046]

[0047] Next, the pre-weighted cosine function of the cone-beam X-ray pair projection data is equivalent to the pre-weighting factor reconstructed under the imaging geometry of the parallel fan beam. Let the angle between the line connecting the X-ray source S' and the center O of the virtual cylindrical detector and the projection point P of the X-ray source S' on the virtual cylindrical detector be α, and the projection of P on the circular orbit plane be P'. First, calculate the value of S'P'.

[0048]

[0049] Next, calculate the value of S'P.

[0050]

[0051] Finally, calculate the cosine function according to the ratio of S'P' to S'P.

[0052]

[0053] Use the cosine function cosα to pre-weight the projection data obtained in step 3, and we have

[0054] p1(s, c, θ) = p0(s, c, θ)cosα (13)

[0055] Step 5: Construct a three-dimensional empirical weighting function based on the conjugate ray method to perform secondary weighting on the pre-weighted data in step 4;

[0056] As Figure 7 shown, A is the reconstruction point on the workpiece to be detected. The projection of the ray emitted by the X-ray source S' passing through A on the circular scanning orbit plane is the straight line l. The intersection point of the straight line l and the circular scanning orbit different from S' is S C ’, S’A and S C ’A are conjugate rays to each other. Denote the angle between S’A and l as α, and denote the angle between S C ’A and l as α C , the distance from P to the circular scanning orbit plane is d, and α C can be considered as a function of α. To make the reconstruction result more reliable, a larger weight is given to the ray with a smaller cone angle in the conjugate rays. This process of increasing the weight first requires defining the difference degree between the conjugate rays S’A and S C ’A

[0057]

[0058] After that, according to the S-shaped curve property of the sigmoid function image, construct an empirical three-dimensional weighting function for α

[0059]

[0060] Then, in order to increase the weight difference between conjugate rays as the distance d from P to the circular orbit plane increases, an empirical constant k greater than 0 is set, and finally a three-dimensional weighting function of S’A is obtained.

[0061]

[0062] Use the three-dimensional weighting function w 3D (α, d) to perform secondary weighting on the projection data obtained in step 4, and there is

[0063] p2(s, c, θ) = p1(s, c, θ)w 3D (α, d) (17)

[0064] Step 6: Use the fast Fourier transform to implement the convolution operation of the one-dimensional ramp filter on the data after secondary weighting in step 5, so as to achieve the filtering effect;

[0065] Use the one-dimensional ramp filter to perform a convolution operation on the projection data obtained in step 4, and there is

[0066] p3(s, c, θ) = p2(s, c, θ) * h(s) (18)

[0067] In the formula, h(s) is the one-dimensional ramp filter, and the filtering of p2(s) is performed along the s direction. The actual filtering path is a curve at the same height on the virtual cylindrical detector. The convolution operation can be implemented by the fast Fourier transform, that is, through the fast Fourier transform, p2(s) is converted into P2(ω), and h(s) is converted into H(ω), and H(ω) satisfies

[0068] H(ω) = |ω| (19)

[0069] After that, P2(ω) and H(ω) are multiplied to obtain P3(ω), and then the inverse fast Fourier transform is performed on P3(ω) to obtain p3(s). In order to eliminate signal aliasing caused by the periodicity of the Fourier transform and reduce artifacts in image reconstruction, zero-padding is used before the convolution operation to double the length of the filter, and the doubling factor is 2 or 4.

[0070] Step 7: Perform a back-projection operation on the data after filtering in step 6 to obtain the value of the reconstruction point;

[0071] Perform a back-projection operation on the data after filtering in step 6 according to the following formula

[0072]

[0073] Wherein, f(x, y, z) is the value of the reconstructed point. For the back-projection operation, considering the property that the operations between various projections are not correlated, cluster technology is used for parallel projection angles, and the projections at different angles are distributed to different nodes in the cluster. After the calculation is completed, the control node accumulates the calculation results of each computing node and then outputs them.

[0074] Step 8: Operate on all the points to be reconstructed according to Steps 4 - 7, and finally obtain the 3D reconstruction result.

[0075] Figure 8 is a physical diagram of a standard spherical plate. After obtaining the projection data by scanning the workpiece to be detected using the method of Step 1, determine the number and position of the reconstructed points, and reconstruct all the reconstructed points according to Steps 4 - 7 to obtain voxel data. During the reconstruction process, the method of constant memory is used to save other intermediate variables that are only related to the projection angle and trigonometric functions. At the same time, each thread of the computer independently completes the reconstruction tasks of a series of points in the corresponding z direction, thereby reducing repeated calculations. After obtaining all the voxel data, trilinear interpolation is used on it to improve the resolution, and finally the 3D reconstruction result of the standard spherical plate is obtained, as Figure 9 shown. Thus, the 3D image reconstruction of industrial CT is achieved, and that's it.

[0076] The above is only a preferred specific implementation manner of the method of the present invention, but the protection scope of the method of the present invention is not limited thereto. Any changes or substitutions that can be easily thought of by those skilled in the art within the technical scope disclosed by the method of the present invention should be covered within the protection scope of the method of the present invention. Therefore, the protection scope of the method of the present invention should be subject to the protection scope of the claims.

Claims

1. An industrial CT three-dimensional image reconstruction method based on data rearrangement and conjugate rays, characterized in that, The implementation is specifically carried out according to the following steps: Step 1: Detect the workpiece using the circular orbit scanning mode for the cone-beam X-ray of the industrial CT to obtain the projection data of the flat detector at different scanning angles; Step 2: Perform a homothetic transformation on the flat detector in Step 1 to obtain a virtual flat detector smaller than the real detector, transfer the projection data on the real detector to the virtual flat detector, then convert the three-dimensional cone-beam X-ray in Step 1 into multiple two-dimensional virtual cone-beam X-ray sources, and further convert the virtual flat detector into a virtual cylindrical detector; Step 3: Perform bilinear interpolation of four neighboring points on the projection data on the virtual flat detector in Step 2, and fill the interpolated projection data into the virtual cylindrical detector in Step 2; Step 4: Perform pre-weighting on the projection data obtained in Step 3 using the cosine function of the X-ray incident angle; Step 5: Construct a three-dimensional empirical weighting function using the method based on conjugate rays, and perform secondary weighting on the pre-weighted data in Step 4. The secondary weighting process is as follows: Let A be the reconstructed point. The projection of S'A on the orbital plane is l, and the intersection of l and the orbit different from S' is S C ’. Denote the angle between S'A and l as α, and denote S C ’A and l as α C . The distance from P to the orbital plane is d. Define the difference degree between S'A and S C ’A After that, construct the function Then, set a constant k greater than 0, and finally obtain the three-dimensional weighting function Finally, use w 3D (α, d) to perform secondary weighting on the pre-weighted projection data p2(s, c, θ) = p1(s, c, θ) w 3D (α, d) P represents the projection on the virtual flat detector; s represents the direction perpendicular to the line connecting the X-ray source and the center of the virtual cylindrical detector; c represents the vertical direction of the virtual flat detector; θ represents the angle of counterclockwise rotation of the s direction relative to the x direction of the scanning object coordinate system; Step 6: Use the fast Fourier transform to implement the convolution operation of the one-dimensional ramp filter and the data after secondary weighting in Step 5, so as to achieve the filtering effect; Step 7: Perform back-projection operation on the filtered data in Step 6 to obtain the value of the reconstructed point; Step 8: Perform operations on all points to be reconstructed according to Steps 4 - 7 to obtain the three-dimensional reconstruction result.

2. The industrial CT three-dimensional image reconstruction method according to claim 1, wherein The process of converting the real detector into a virtual detector in Step 2 is as follows: Taking the X-ray source as the homothetic center, perform a homothetic transformation on the real flat detector to make it coincide with the center plane of the rotary table to obtain the virtual flat detector, and transfer the projection data on the real detector to the virtual flat detector Then, establish the conversion relationship between the virtual flat detector and the virtual cylindrical detector x represents the x direction of the scanning object coordinate system; y represents the y direction of the scanning object coordinate system; a represents the horizontal direction of the virtual flat detector; b represents the direction of the line connecting the X-ray source and the center of the virtual flat detector; Indicates the angle by which the virtual plane detector rotates counterclockwise relative to the rotating table; R represents the distance from the X-ray source to the center of the virtual flat detector.

3. The industrial CT three-dimensional image reconstruction method according to claim 1, characterized in that The filling process of the virtual cylindrical detector in Step 3 is as follows: Let the projection of a point on the virtual plane detector be The projection of the adjacent point diagonally down and to the left of this point is The projection of the adjacent point diagonally up and to the left of this point is The projection of the adjacent point diagonally down and to the right of this point is The projection of the adjacent point diagonally up and to the right of this point is Calculate the difference between p and p 11 ​ u = a - a0 v = c - c0 For Using four-neighbor bilinear interpolation, there is Use to fill the virtual cylindrical detector, and make the number of filling points at each height be an integer power of 2. The filling relationship is 4. The industrial CT three-dimensional image reconstruction method according to claim 1, wherein The cosine weighting process in Step 4 is as follows: Calculate the height of the projection point on the virtual cylindrical detector from the triangle similarity relationship Let the angle between the line connecting the X-ray source S’ and the center of the virtual cylindrical detector and the line from S’ to the projection point P on the virtual cylindrical detector be α, and the projection of P on the circular orbit plane is P’. First calculate S’P’, then calculate S’P, and finally calculate the cosine function Use cosα to perform pre-weighting on the interpolated projection data, and there is p1(s, c, θ) = p0(s, c, θ)cosα h represents the height of the projection point on the virtual cylindrical detector; t represents the direction of the line connecting the X-ray source S' and the center of the virtual cylindrical detector; z represents the z direction of the coordinate system of the scanned object.

5. The industrial CT three-dimensional image reconstruction method according to claim 1, characterized in that The filtering process in step 6 is as follows: Perform a convolution operation on the projection data after quadratic weighting using a one-dimensional ramp filter h(s), where p3(s, c, θ) = P2(s, c, θ) * h(s) Before the convolution operation, use zero-padding to double the length of the filter. The doubling factor is 2 or 4, and the convolution operation is implemented through the fast Fourier transform and the inverse fast Fourier transform.

6. The industrial CT three-dimensional image reconstruction method according to claim 1, wherein The backprojection process in step 7 is as follows: Perform a backprojection operation on the filtered data. Considering the property that the operations between projections are not correlated, use the cluster technology for parallel projection angles, distribute the projections at different angles to different nodes in the cluster, and after the calculation is completed, the control node accumulates the calculation results of each computing node and outputs them.

Images