Imaging geometry online calibration method

Abstract

The invention discloses an imaging geometry online calibration method, which comprises the following steps that: a marker motif is scanned by a high-degree-of-freedom CT (Computed Tomography) to obtain a first group of projection drawings, and the marker motif comprises an upper half part and a bottom part; processing the projection data based on the first group of projection drawings to obtain initial geometric information; removing the upper half part of the marker motif, and scanning the bottom part of the marker motif and the measured object by the high-degree-of-freedom CT to obtain a second group of projection drawings; correcting the initial geometric information based on a global optimization algorithm and the second group of projection drawings to obtain final optimized geometric information; three-dimensional reconstruction is carried out based on final optimized geometric information, the marker motif device used in the method occupies little imaging view and is located at the bottom of the view, imaging of a measured object is not interfered, and the engineering applicability is high.

Description

Technical Field

[0001] The invention belongs to the technical field of X-ray computed tomography (CT) imaging, and in particular relates to an imaging geometry online calibration method. Background Art

[0002] X-ray-based computed tomography (CT) is one of the current mainstream non-destructive testing technologies. Figure 1 As shown in the figure, the basic principle is that X-rays penetrate the object being tested. Due to the different absorption rates of different materials for radiation, the X-rays that pass through the object will leave counts of different intensities on the pixels at different positions of the detector. All the detector pixel counts are arranged together to form a projection image. CT obtains projection images at different angles by rotating the radiation source and the detector relative to the object being tested on the stage, and then uses reconstruction algorithms (such as FDK) to reconstruct the internal information of the object.

[0003] As a new type of CT structure, high-degree-of-freedom CT mainly uses a flexible motion device to load the imaging components X-ray source and flat-panel detector. During imaging, the driving device rotates the ray source and detector around the object to be measured to obtain projection data at different angles, thereby reconstructing the tomographic profile of the object to be measured, and non-destructively observing the fine structure inside the object to be measured. The characteristics of the high-degree-of-freedom CT system are: 1. Strong flexibility, high degree of freedom, and multiple imaging modes; 2. It can not only perform large-field imaging of large-sized objects, but also achieve high-resolution imaging of local positions of complex structure objects, which is suitable for scanning more types of objects.

[0004] Due to the flexibility of its structure, high-degree-of-freedom CT may have non-circular acquisition trajectories, asynchronous motion of the detector and the radiation source, and low precision of the motion of each axis, which may affect the quality of the reconstructed image, causing quality degradation, distortion or artifacts, and failure to achieve the purpose of detection. Accurate and effective imaging geometry calibration is one of the key technologies for high-degree-of-freedom CT imaging, which can ensure the accuracy of data acquisition.

[0005] Existing calibration methods are mainly based on offline correction methods of marked phantoms, such as Figure 2 As shown, a standard marker phantom (such as a marker ball) is placed in the field of view for a CT scan, and the projection geometric parameters at a single angle are expressed as a homogeneous matrix using a correction method based on the projection matrix of the correction phantom. The projection matrix at the corresponding angle is calculated through the one-to-one correspondence between the marker points in the phantom at a single angle and the corresponding two-dimensional projection points. In the process of back-projection reconstruction, the corresponding back-projection position information is calculated using the geometric information contained in the projection matrix, so that the position information is used for reconstruction of the actual object being measured.

[0006] The above-mentioned offline calibration method is applicable to the case where the high-degree-of-freedom CT acquisition system has a very high repetitive positioning accuracy of the motion axis. However, in reality, when the motion axis of the high-degree-of-freedom CT acquisition system repeatedly moves, the geometric information of the two repeated movements cannot be completely consistent, and the deviation between the two depends on the repetitive positioning accuracy of the motion axis of the high-degree-of-freedom CT acquisition system. When the accuracy cannot meet the requirements of high-resolution imaging accuracy, that is, when the error of the positioning accuracy is greatly different from the size of the reconstructed pixel, it may not be possible to use the geometric information obtained from a single marker scan to reconstruct the projection data of each subsequent scan. At this time, the offline calibration method is no longer applicable, and it is necessary to study an online correction method that can be applied to real-time acquisition. Summary of the invention

[0007] In order to solve the above technical problems, the present invention provides an imaging geometry online calibration method, which for the first time proposes a feasible solution to the problem of high-degree-of-freedom CT imaging geometry online calibration, and obtains good correction effect through algorithm verification.

[0008] To achieve the above object, the present invention adopts the following technical solutions: A method for online calibration of imaging geometry, the method comprising: Step 1: Scanning the marker phantom with a high degree of freedom CT to obtain a first set of projection images, wherein the marker phantom includes an upper part and a bottom part, and both the upper part and the bottom part include a standard small ball made of a metal material; Step 2: Processing the projection data based on the first set of projection images to obtain initial geometric information; Step 3, removing the upper part of the marker phantom, and scanning the bottom part of the marker phantom and the object to be measured by the high degree of freedom CT to obtain a second set of projection images; Step 4: Based on the global optimization algorithm and the second set of projection images, the initial geometric information is corrected to obtain the final optimized geometric information; Step 5: Perform 3D reconstruction based on the final optimized geometric information.

[0009] Furthermore, the upper part and the bottom part of the marker model are mounted together by a rotating buckle, and the standard balls are distributed at equal intervals along a spiral line.

[0010] Furthermore, the step 2 includes: identifying the center of the standard sphere on the first group of projection images, and obtaining the position coordinates pt of the center of each standard sphere in the projection image in the detector coordinate system: , It represents the coordinates of the center of the standard ball on one of the projections in the first group of projections; Define the geometric information calculated for each projection as twelve-dimensional: , Among them, S, D, U, V are three-dimensional geometric parameters defined in the three-dimensional reconstruction coordinate system, U and V are two vectors that determine the detector posture and whose modulus is the length of one detector pixel, S is the coordinate of the ray source, and D is the coordinate of the lower left pixel of the detector; Define vpos as the three-dimensional coordinates of the center of the standard ball in the marker phantom in the three-dimensional reconstruction coordinate system: , The projection matrix P of the center of the standard ball projected onto the detector coordinate system is defined as: , , In the formula, the superscript T represents the transpose of the matrix; Set the intermediate vector w to: , Calculate the projection point coordinates pt_proj of the center of the standard ball in the detector coordinate system: , Based on the assumption that the projection point coordinate pt_proj and the position coordinate pt are the same point, initial geometric information is obtained by calculating each projection image of the first group of projection images.

[0011] Furthermore, the step 3 includes: separating the upper part of the marker phantom and moving it out of the field of view, placing the object to be measured on the bottom part of the marker phantom, and performing a second scan with the high degree of freedom CT to obtain a second set of projection images.

[0012] Further, the step 4 comprises: Obtain the center position pt' of each standard ball of the bottom part of the phantom in the second group of projections projected onto the detector coordinate system; The geometric information of the second group of projections that are generated in the same time order as the first group of projections is used as input, and the data processing method of step 2 is used again to calculate the projection point coordinates pt_proj' of the center of the standard ball of the bottom part of the phantom on the detector coordinate system, and calculate the Euclidean distance between pt' and pt_proj' : , The target function ObjectFunc is obtained by summing the Euclidean distances pterror calculated for all standard balls in the bottom part: , For each projection image in the second group of projection images, the initial value of the geometric information dimension with poor repeatability of the high-degree-of-freedom CT system is selected as the input variable of the optimization algorithm, and the variation range of each variable is specified. The initial values ​​of other geometric information dimensions are used as input constants. The values ​​of the variables are changed in each iteration, and the objective function value is calculated. When the objective function value converges to below a certain value, the changed geometric information is output to obtain the final optimized geometric information.

[0013] Furthermore, step 5 includes performing three-dimensional reconstruction on the second group of projection images using the FDK algorithm according to the final optimized geometric information and the initial geometric information.

[0014] The beneficial effects of the present invention are: The present invention proposes an online calibration method for imaging geometry, which is based on global iterative optimization of key parameters of high-degree-of-freedom CT. The marker phantom device used occupies very little imaging field of view and is located at the bottom of the field of view, without any interference with the imaging of the object being measured. The method is suitable for online calibration of high-precision imaging of CT systems. BRIEF DESCRIPTION OF THE DRAWINGS

[0015] Figure 1 It is a schematic diagram of the structure of a CT system in the prior art; Figure 2 It is a schematic diagram of the principle of the marker-based CT geometric information analysis method in the prior art; Figure 3 It is a schematic diagram of the standard marker motif; Figure 4 This is a flow chart of a CT imaging geometry online calibration method of the present invention; Figure 5 It is a schematic diagram of the phantom projection image; Figure 6 is a schematic diagram of the detector coordinate system; Figure 7 Schematic diagram of the detector's deflection in different directions; Figure 8 A schematic diagram showing the comparison of reconstruction results. DETAILED DESCRIPTION

[0016] The present invention will be further described below in conjunction with the accompanying drawings and embodiments.

[0017] The present invention provides a marker phantom specific structure used in an imaging geometry online calibration method. Figure 3 As shown, the marker phantom includes an upper part and a bottom part, both of which are cylinders containing standard small balls. The two cylinders are mounted together by a rotating buckle. The cylinder is generally made of resin material, and the standard small ball is made of metal material, such as Figure 3As shown, the upper cylinder contains 7-9 metal balls, and the bottom cylinder contains 2-4 metal balls, wherein the distribution requirement of the metal balls is: any two circles in the projection image obtained by CT scanning shall not block each other, and the specific arrangement depends on the CT field of view, and the general form is that the metal balls are spirally equally spaced and ascending. Without loss of generality, the shape of the external structure for loading metal balls is not limited to a cylinder, and the external structure does not affect the specific implementation of the invention. Before imaging, the marker phantom needs to be calibrated in three coordinates, and a three-dimensional coordinate system of the marker phantom is established to obtain the spherical center coordinates of all standard spheres in the marker phantom. The calibration process is: during imaging, the marker phantom is scanned by high-degree-of-freedom CT, and then the spherical center coordinates of all calibrated standard spheres and the circle center coordinates in the projection image are combined to analyze and obtain the geometric information of the imaging. During online calibration, keep the bottom part of the phantom still, remove the upper part, place the object to be measured on the bottom part, and use the high-degree-of-freedom CT to scan the bottom part and the object to be measured at the same time. The center coordinates of the standard ball in the bottom part and the center coordinates of the circle in the projection are used as input data, and online calibration is performed through the corresponding algorithm. The advantages of this device are: before calibration, the entire marker phantom is scanned only once, which does not bring about cumbersome calibration steps; during online real-time calibration, only the bottom part is involved in imaging, which occupies very little field of view, and is at the bottom of the field of view, so it will not interfere with the imaging of the object to be measured, and it has stronger engineering applicability.

[0018] like Figure 4 As shown, the present invention provides a schematic diagram of an imaging geometry online calibration method. First, the entire marker phantom is scanned and analyzed to obtain a set of initial geometric information. The geometric information includes: the coordinates of a certain pixel of the detector in the three-dimensional coordinate system of the marker phantom calibration, the coordinates of the ray source, and two vertical vectors that determine the surface of the detector pixel array. During online calibration, the upper half of the phantom is separated, and the object to be measured can be placed on the bottom part. On the basis of the previous initial geometric information, the projection of this scan and the coordinates of the center of the standard sphere in the bottom part are used as input to set the optimized objective function value, and the initial geometric information is modified by the corresponding algorithm to obtain imaging geometric information that is more suitable for this scan. The algorithm used is a global optimization algorithm, which aims to find the global optimal value of the function, and is usually used to deal with optimization problems with multiple local minima. Specifically: Step 1: Scan the entire marker phantom with high-degree-of-freedom CT; The geometric system of the high-degree-of-freedom CT imaging system can flexibly and accurately image the object being measured. The ray source and detector can complete the projection data acquisition by rotating around the object being measured, and then the 3D data of the object being measured can be obtained by using the CT reconstruction algorithm. Under this CT system, the above-mentioned marker phantom can be scanned to obtain the first set of projection images.

[0019] Step 2: Processing the projection data based on the first group of projection images; The center of the standard ball is identified on the first set of projections, where the projections are as follows: Figure 5 As shown. Get the position coordinates pt of the center of each standard ball in the projection image in the detector coordinate system: , It represents the coordinates of the center of the standard ball on one of the projections in the first group of projections; The detector coordinate system diagram is as follows: Figure 6 As shown, each point represents a pixel in the detector.

[0020] The geometric information calculated for each projection is defined as twelve-dimensional, specifically: , Among them, S, D, U, V are three-dimensional geometric parameters defined in the three-dimensional reconstruction coordinate system, U, V are two vectors that determine the detector posture and whose modulus is the length of a detector pixel (expressed as u, v in the detector coordinate system, such as Figure 7 As shown in the figure), S is the coordinate of the ray source, and D is the coordinate of the lower left pixel of the detector.

[0021] Define vpos as the three-dimensional coordinates of the center of the standard ball in the marker phantom in the three-dimensional reconstruction coordinate system: , The projection matrix P of the sphere center onto the detector coordinate system is defined as: , , Where the superscript T represents the transpose of the matrix.

[0022] Set the intermediate vector w to: , Calculate the projection point coordinates pt_proj of the sphere center in the detector coordinate system: , For seven or more spheres in the model, an analytical formula is established, assuming that the above projection point coordinates pt_proj and pt are the same point, that is, the geometric information (S, D, U, V) corresponding to the projection image can be calculated through the projection matrix relationship between the center of the projection image and the center of the sphere. Calculate all the projection images in step 1 to obtain the initial geometric information.

[0023] Step 3: High-degree-of-freedom CT scans the bottom part of the marker phantom and the object to be measured; The upper part of the marker phantom is separated and moved out of the field of view, while the bottom part of the phantom remains stationary. The object to be tested is placed on the bottom part of the marker phantom, and the CT scan is performed again for the second time to obtain a second set of projection images.

[0024] Step 4: Correcting the initial geometric information based on a global optimization algorithm; Select the dimension of geometric information that needs to be changed: Select the dimension that may change from the initial geometric information in this scan and has a significant impact on the quality of the reconstructed image from the twelve-dimensional geometric information. The selection criteria are the dimensions with more obvious statistical changes through repeated calculation of geometric information experiments on the CT system. The dimensions of geometric information with poor repeatability may be different for each CT system, depending on the specific equipment. The changes in the two posture vectors U and V of the detector are abstracted into changes in three posture angles. The posture angles are as follows: Figure 7 As shown in the figure, the first attitude angle θ is the angle at which the U vector rotates around the V vector; the second attitude angle : The angle at which the V vector rotates around the U vector; the third attitude angle :The angle at which the U and V vectors rotate around the normal vector of the face they are on.

[0025] Get the position pt' of each center of the bottom part of the phantom in the detector coordinate system in the second group of projections, Take the geometric information of the second set of projections and the first set of projections that have the same generation order as the first set of projections as input, and process them again according to the data processing method in step 2 to calculate the projection point pt_proj' of the sphere center of the bottom part of the model on the detector coordinate system. Calculate the Euclidean distance between pt' and pt_proj' : , The objective function ObjectFunc is obtained by summing the pterrors calculated for the n standard balls at the bottom. The ultimate goal of the optimization is to make the value of the objective function closer to 0: , The optimization algorithm is a global optimization algorithm. The initial values ​​of the geometric information dimensions with poor repeatability of the CT system are used as input variables of the optimization algorithm, and the range of variation of each variable is specified. The initial values ​​of other dimensions are used as input constants. The values ​​of the variables are changed in each iteration, and the objective function value is calculated. During the iteration, when the objective function value converges to below a certain value, the changed geometric information is output. This operation is repeated for each projection in the second group of projections to obtain the final optimized geometric information.

[0026] Step 5: Perform 3D reconstruction based on the final optimized geometric information; The second set of projection images are reconstructed using the FDK algorithm based on the final optimized geometric information and the initial geometric information. The reconstruction results before and after correction are compared. Figure 8 As shown, the left picture is the reconstruction of the initial geometric information before correction, and the right picture is the reconstruction of the geometric information after correction. It can be seen that the edge of the reconstructed image with the corrected geometric information is clearer.

[0027] The specific embodiments described above further illustrate the objectives, technical solutions and beneficial effects of the present invention in detail. It should be understood that the above description is only a specific embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principles of the present invention should be included in the scope of protection of the present invention.

Claims

1. An imaging geometry online calibration method, characterized in that: The method comprises: Step 1: Scanning the marker phantom with a high degree of freedom CT to obtain a first set of projection images, wherein the marker phantom includes an upper part and a bottom part, and both the upper part and the bottom part include a standard small ball made of a metal material; Step 2: Processing the projection data based on the first set of projection images to obtain initial geometric information; Step 3, removing the upper part of the marker phantom, and scanning the bottom part of the marker phantom and the object to be measured by the high degree of freedom CT to obtain a second set of projection images; Step 4: Based on the global optimization algorithm and the second set of projection images, the initial geometric information is corrected to obtain the final optimized geometric information; Step 5: Perform 3D reconstruction based on the final optimized geometric information.

2. The method for online calibration of imaging geometry according to claim 1, characterized in that: In the step 1, the upper part and the bottom part of the marker model are mounted together by a rotating buckle, and the standard balls are distributed at equal intervals along a spiral line.

3. The method for online calibration of imaging geometry according to claim 1, characterized in that: The step 2 includes: identifying the center of the standard sphere on the first group of projection images, and obtaining the position coordinates pt of the center of each standard sphere in the projection image in the detector coordinate system: ， It represents the coordinates of the center of the standard ball on one of the projections in the first group of projections; Define the geometric information calculated for each projection as twelve-dimensional: ， Among them, S, D, U, V are three-dimensional geometric parameters defined in the three-dimensional reconstruction coordinate system, U and V are two vectors that determine the detector posture and whose modulus is the length of one detector pixel, S is the coordinate of the ray source, and D is the coordinate of the lower left pixel of the detector; Define vpos as the three-dimensional coordinates of the center of the standard ball in the marker phantom in the three-dimensional reconstruction coordinate system: ， The projection matrix P of the center of the standard ball projected onto the detector coordinate system is defined as: ， ， In the formula, the superscript T represents the transpose of the matrix; Set the intermediate vector w to: ， Calculate the projection point coordinates pt_proj of the center of the standard ball in the detector coordinate system: ， Based on the assumption that the projection point coordinate pt_proj and the position coordinate pt are the same point, initial geometric information is obtained by calculating each projection image of the first group of projection images.

4. The method for online calibration of imaging geometry according to claim 1, characterized in that: The step 3 comprises: separating the upper part of the marker phantom and moving it out of the field of view, placing the object to be measured on the bottom part of the marker phantom, and performing a second scan on the high-degree-of-freedom CT to obtain a second set of projection images.

5. The method for online calibration of imaging geometry according to claim 1, characterized in that: The step 4 comprises: Obtain the center position pt' of each standard ball of the bottom part of the phantom in the second group of projections projected onto the detector coordinate system; The geometric information of the second group of projections that are generated in the same time order as the first group of projections is used as input, and the data processing method of step 2 is used again to calculate the projection point coordinates pt_proj' of the center of the standard ball of the bottom part of the phantom on the detector coordinate system, and calculate the Euclidean distance between pt' and pt_proj' : ， The target function ObjectFunc is obtained by summing the Euclidean distances pterror calculated for all standard balls in the bottom part: ， For each projection image in the second group of projection images, the initial value of the geometric information dimension with poor repeatability of the high-degree-of-freedom CT system is selected as the input variable of the optimization algorithm, and the variation range of each variable is specified. The initial values ​​of other geometric information dimensions are used as input constants. The values ​​of the variables are changed in each iteration, and the objective function value is calculated. When the objective function value converges to below a certain value, the changed geometric information is output to obtain the final optimized geometric information.

6. The method for online calibration of imaging geometry according to claim 1, characterized in that: The step 5 includes performing three-dimensional reconstruction on the second group of projection images using the FDK (Feldkamp, ​​Davis and Kress) algorithm according to the final optimized geometric information and the initial geometric information.

Images