Method and device for calibrating geometry parameters of CT system based on reconstruction feedback optimization

Abstract

The application provides a CT system geometry parameter calibration method and device based on reconstruction feedback optimization, and the method comprises the following steps: establishing a CT system geometry coordinate system and determining geometry parameters to be calibrated; setting a single-layer ball array calibration board; controlling the CT system to collect a projection image set; iteratively performing the following operations until a stop condition is met: taking the position of each particle in a particle swarm as a set of candidate geometry parameters, performing three-dimensional reconstruction according to the projection image set and the corresponding candidate geometry parameters under the constraint of the CT system geometry coordinate system, and obtaining three-dimensional reconstruction volume data corresponding to each set of candidate geometry parameters; extracting sphere data from the three-dimensional reconstruction volume data, solving a target function value related to an error according to the sphere data, and updating the position of the particle; and taking the global optimal solution obtained when the stop condition is met as the final calibrated geometry parameter. The method can improve the accuracy of CT geometry calibration.

Description

Technical Field

[0001] This invention relates to computed tomography (CT) technology, and more particularly to a method and apparatus for calibrating geometric parameters of a CT system based on reconstruction feedback optimization. Background Technology

[0002] In the field of industrial nondestructive testing, as product structures become increasingly complex, the scope of defect detection has expanded from surface defects to internal structural defects. For products such as castings, weldments, composite materials, electronic packaging components, lithium batteries, and precision assemblies, defects such as porosity, inclusions, cracks, delamination, incomplete soldering, internal fractures, and assembly misalignments are often hidden within the object, making them difficult to reliably identify using conventional two-dimensional visual inspection. Therefore, X-ray CT technology, capable of acquiring three-dimensional structural information of the object's internal structure, has become an important tool in industrial quality inspection and nondestructive evaluation.

[0003] Computed tomography (CT) technology acquires three-dimensional volumetric data of an object by acquiring X-ray projection data from multiple angles and using reconstruction algorithms to reconstruct the attenuation coefficient distribution within the object. The results of CT reconstruction depend not only on the X-ray source, detector, number of projections, and the reconstruction algorithm itself, but also heavily on the accuracy of the system's geometric model. The reconstruction algorithm needs to accurately describe the spatial relationships between the X-ray source, rotation axis, detector plane, and detector pixel coordinates at each projection angle. If the actual system geometry differs from the geometric parameters used in the reconstruction algorithm, the path calculation of the rays in voxel space will deviate, leading to structural ghosting, blurred edges, dimensional distortion, deformation of circular targets, and defects in the reconstructed object.

[0004] Patent CN112634353B proposes a self-calibration method for a geometric calibration phantom in a CBCT system. The core idea is to acquire the calibration phantom and determine its approximate spatial coordinates, collect projection images of the phantom in various orientations, determine the coordinates of the marker's projection center, calculate the system's geometric parameters for each orientation, and construct constraints and an objective function based on the projection center coordinates and system geometric parameters to calculate the marker's target spatial coordinates. This scheme primarily addresses the problem of self-calibrating the marker's spatial coordinates when there are processing errors in the marker's actual spatial position within the calibration phantom. This method mainly utilizes the relationship between the marker's center position in two-dimensional projection and geometric parameters for solution. However, two-dimensional projection errors are not always completely consistent with the final three-dimensional reconstruction quality. For example, certain combinations of geometric parameters may show similar two-dimensional projection errors, but in the three-dimensional reconstruction results, they may produce varying degrees of spherical deformation, edge diffusion, spatial scale distortion, or interlayer misalignment. When geometric parameter errors mainly manifest as deformation, ghosting, spherical distortion, or local structural blurring in the three-dimensional reconstruction, this scheme cannot be directly optimized based on the three-dimensional reconstruction quality. Summary of the Invention

[0005] This invention provides a method and apparatus for calibrating the geometric parameters of a CT system based on reconstruction feedback optimization, which improves the quality of three-dimensional reconstructions by optimizing the geometric parameters of the CT system.

[0006] A method for calibrating geometric parameters of a CT system based on reconstruction feedback optimization, comprising: Establish the geometric coordinate system of the CT system and determine the geometric parameters to be calibrated; Set up a single-layer ball array calibration board; The CT system is controlled to acquire projection images of a single-slice spherical array calibration plate at different projection angles according to the scanning trajectory, forming a set of projection images; The following operations are performed iteratively until the stopping condition is met: the position of each particle in the particle swarm is taken as a set of candidate geometric parameters. Under the constraints of the geometric coordinate system of the CT system, three-dimensional reconstruction is performed based on the set of projected images and the corresponding candidate geometric parameters to obtain the three-dimensional reconstructed volume data corresponding to each set of candidate geometric parameters; sphere data is extracted based on the three-dimensional reconstructed volume data, and the objective function value with respect to the error is solved based on the sphere data to update the position of the particles; The globally optimal solution obtained when the stopping condition is met is used as the final calibration geometric parameters.

[0007] Furthermore, in the geometric coordinate system of the CT system, the origin is the intersection of the plane containing the rotation axis and the center of rotation, the Z-axis is the rotation axis, and the X-axis and Y-axis are located in a plane perpendicular to the rotation axis and form a right-handed coordinate system with the Z-axis; In the geometric coordinate system of the CT system, the geometric parameters to be calibrated include the three-dimensional position coordinates of the X-ray source under the initial projection angle, the three-dimensional position coordinates of the detector center under the initial projection angle, the three-dimensional displacement vector per unit pixel in the X direction of the detector under the initial projection angle, and the three-dimensional displacement vector per unit pixel in the Y direction of the detector under the initial projection angle.

[0008] Furthermore, the single-layer sphere array calibration plate includes a substrate and spheres arranged in an array on the substrate.

[0009] Furthermore, 3D reconstruction is performed based on the set of projected images and the corresponding candidate geometric parameters, including: Based on the candidate geometric parameters and the projection angle of each projection image in the projection image set, generate equivalent geometric parameters for the projection angle of each projection image. The three-dimensional coordinates of the pixels in the projected image corresponding to the detector space are calculated based on the equivalent geometric parameters of the projection angle. The normal vector of the detector plane and the normal distance from the ray source to the detector plane are calculated based on the equivalent geometric parameters of the projection angle. The corresponding projected image is weighted by cone-beam averaging based on the pixel's three-dimensional coordinates and the normal distance from the X-ray source to the detector plane to obtain a weighted projected image. The weighted projection image is filtered along the X direction of the detector to obtain a filtered weighted projection image; Based on the normal vector under the projection angle and the normal distance from the X-ray source to the detector plane, the filtered weighted projection image is accumulated by voxel projection and back projection to obtain the three-dimensional reconstructed volume data corresponding to the candidate geometric parameters.

[0010] Furthermore, the three-dimensional position coordinates of the ray source, the three-dimensional position coordinates of the detector center, the three-dimensional displacement vector of the unit pixel in the X direction of the detector, and the three-dimensional displacement vector of the unit pixel in the Y direction of the detector in the current candidate geometric parameters are multiplied by the rotation matrix corresponding to the projection angle to obtain the equivalent geometric parameters of the projection angle. The equivalent geometric parameters of the projection angle include the three-dimensional position coordinates of the X-ray source under the projection angle, the three-dimensional position coordinates of the detector center under the projection angle, the three-dimensional displacement vector per unit pixel in the X direction of the detector under the projection angle, and the three-dimensional displacement vector per unit pixel in the Y direction of the detector under the projection angle.

[0011] Further, calculating the three-dimensional coordinates of pixels in the projected image corresponding to the detector space based on the equivalent geometric parameters of the projection angle includes: Calculate the center pixel coordinates of the detector based on its width and height; Calculate the offset of each pixel in the projected image relative to the coordinates of the detector's center pixel; Based on the offset, the three-dimensional displacement vector of a unit pixel in the X direction of the detector under the projection angle, and the three-dimensional displacement vector of a unit pixel in the Y direction of the detector under the projection angle, the displacement vector of the pixel in the detection space is calculated. Based on the three-dimensional position coordinates of the detector center under the projection angle and the displacement vector corresponding to the pixel in the detection space, the three-dimensional coordinates of each pixel in the projected image in the detector space are obtained.

[0012] Further, calculating the normal vector under the projection angle and the normal distance from the ray source to the detector plane based on the equivalent geometric parameters of the projection angle includes: The normal vector under the projection angle is calculated based on the three-dimensional displacement vector of a unit pixel in the X direction of the detector under the projection angle and the three-dimensional displacement vector of a unit pixel in the Y direction of the detector under the projection angle. Based on the normal vector, the three-dimensional position coordinates of the X-ray source under the projection angle, and the three-dimensional position coordinates of the detector center under the projection angle, the normal distance from the X-ray source to the detector plane is calculated.

[0013] Furthermore, cone-beam weighting is applied to the corresponding projected image based on the pixel's three-dimensional coordinates and the normal distance from the ray source to the detector plane, including: The cone-beam weighting coefficients of the corresponding pixels are calculated based on the normal distance from the X-ray source to the detector plane and the three-dimensional coordinates of the pixels. The cone-beam weighting coefficients are applied to the corresponding pixels in the projected image to obtain a weighted projected image.

[0014] Furthermore, based on the normal vector at the projection angle and the normal distance from the ray source to the detector plane, the filtered weighted projection image is accumulated through voxel projection and back projection, including: Based on the normal vector, the three-dimensional position coordinates of the X-ray source under the projection angle, and the pre-determined center coordinates of the voxel, calculate the distance of the voxel relative to the X-ray source in the detector normal direction. If the distance of the voxel relative to the X-ray source in the detector normal direction is greater than 0, then the intersection point of the detector plane corresponding to the voxel is calculated based on the normal distance from the X-ray source to the detector plane, the distance of the voxel relative to the X-ray source in the detector normal direction, the three-dimensional position coordinates of the X-ray source under the projection angle, and the center coordinates of the voxel. The intersection points of the detector planes corresponding to the voxels are converted into continuous pixel coordinates on the filtered weighted projection image; Bilinear interpolation is performed at the continuous pixel coordinates on the filtered weighted projection image to obtain the back-projection image; The back-projected images under all projection angles are weighted and accumulated to obtain the 3D reconstructed volume data corresponding to the candidate geometric parameters.

[0015] Furthermore, sphere data is extracted from the 3D reconstructed volume data, and the objective function value with respect to the error is solved based on the sphere data, including: Based on the three-dimensional reconstructed data, threshold segmentation and three-dimensional connected component analysis are performed to calculate the reconstructed sphere center coordinates and reconstructed sphere radius of the sphere in the single-layer sphere array calibration plate. The sphericity error is calculated based on the reconstructed sphere center coordinates and reconstructed sphere radius, and the radius error is calculated based on the reconstructed sphere radius. Based on the reconstructed sphere center coordinates, establish the sphere array coordinate direction matrix and calculate the sphere center array coordinates; Calculate the standard sphere center coordinates based on the standard sphere center coordinates in the physical coordinate system of the single-layer sphere array calibration plate; The array coordinate error is calculated based on the array coordinates of the sphere center and the standard sphere center-centered coordinates, and the sphere spacing error is calculated based on the standard sphere center-centered coordinates and the reconstructed sphere center coordinates. Define the ball count penalty and parameter regularization terms; The objective function value is obtained by weighted summing of the sphericity error, radius error, array coordinate error, sphere spacing error, sphere quantity penalty term, and parameter regularization term.

[0016] Furthermore, based on the three-dimensional reconstructed data, threshold segmentation and three-dimensional connected component analysis are performed to calculate the coordinates of the reconstructed sphere center and the radius of the reconstructed sphere in the single-layer sphere array calibration plate, including: The three-dimensional reconstructed volume data is segmented according to the set segmentation threshold to obtain binary volume data representing the spherical region; Based on the binary data, a three-dimensional connected component analysis was performed to obtain multiple spherical regions; Based on the voxel center coordinates of each voxel in the sphere region and the corresponding 3D reconstructed volume data, the coordinates of the reconstructed sphere center are calculated by gray-scale weighting. The radius of the reconstructed sphere is obtained by performing least-squares fitting on the boundary points of the sphere region.

[0017] Further, establishing the sphere array coordinate direction matrix based on the reconstructed sphere center coordinates includes: Obtain the coordinates of all reconstructed sphere centers to form a set of sphere centers; Calculate the sphere center covariance matrix based on the sphere center set, perform eigenvalue decomposition on the sphere center covariance matrix, and take the unit eigenvectors corresponding to the two largest eigenvalues ​​as the row and column directions of the sphere array coordinate system; The height direction is obtained by cross-multiplying the unit eigenvectors of the row and column directions. The row, column, and height directions together form the coordinate direction matrix of the spherical array.

[0018] Furthermore, the sphericity error is the sum of squares of the normalized spherical residuals of the sphere, and the normalized spherical residual of a single sphere is the ratio of the root mean square residual of the distance from the boundary point of the sphere region to the center of the reconstructed sphere relative to the radius of the reconstructed sphere to the radius of the reconstructed sphere. The radius error is the mean square error between the reconstructed sphere radius and the standard sphere radius of all spheres; The sphere center array coordinates are the product of the transpose of the sphere array coordinate direction matrix and the reconstructed sphere center coordinates and the difference between the reconstructed sphere center coordinates and the mean of the reconstructed sphere center coordinates; The standard sphere center-centered coordinates are the difference between the standard sphere center coordinates and the mean of the standard sphere center coordinates in the physical coordinate system of the single-layer sphere array calibration plate. The array coordinate error is the mean square error of the difference between the array coordinates of the center of all spheres and the standard sphere center-centered coordinates; the sphere spacing error is the mean square error of the difference between the reconstructed sphere center spacing of all two adjacent spheres and the standard sphere center spacing. The ball count penalty term is the absolute value of the difference between the number of balls detected and the actual number of balls under the current candidate geometric parameters; the parameter regularization term is the sum of squared deviations of the changes between the initial geometric parameters and the current candidate parameters.

[0019] A CT system geometric parameter calibration device based on reconstruction feedback optimization applied to the above method includes: The coordinate establishment module is used to establish the geometric coordinate system of the CT system and determine the geometric parameters to be calibrated. The image acquisition module is used to control the CT system to acquire projected images of the single-layer spherical array calibration plate at different projection angles according to the scanning trajectory, forming a set of projected images; The optimization module iteratively performs the following operations until a stopping condition is met: taking the position of each particle in the particle swarm as a set of candidate geometric parameters, and under the constraints of the geometric coordinate system of the CT system, performing 3D reconstruction based on the set of projected images and the corresponding candidate geometric parameters to obtain 3D reconstructed volume data corresponding to each set of candidate geometric parameters; extracting sphere data based on the 3D reconstructed volume data, solving the objective function value with respect to the error based on the sphere data, and updating the particle positions; The parameter determination module is used to take the globally optimal solution obtained when the stopping condition is met as the final calibration geometric parameters.

[0020] Furthermore, in the geometric coordinate system of the CT system, the origin is the intersection of the plane containing the rotation axis and the center of rotation, the Z-axis is the rotation axis, and the X-axis and Y-axis are located in a plane perpendicular to the rotation axis and form a right-handed coordinate system with the Z-axis; In the geometric coordinate system of the CT system, the geometric parameters to be calibrated include the three-dimensional position coordinates of the X-ray source under the initial projection angle, the three-dimensional position coordinates of the detector center under the initial projection angle, the three-dimensional displacement vector per unit pixel in the X direction of the detector under the initial projection angle, and the three-dimensional displacement vector per unit pixel in the Y direction of the detector under the initial projection angle.

[0021] Furthermore, the single-layer sphere array calibration plate includes a substrate and spheres arranged in an array on the substrate.

[0022] Furthermore, the optimization module performs 3D reconstruction based on the set of projected images and the corresponding candidate geometric parameters, including: Based on the candidate geometric parameters and the projection angle of each projection image in the projection image set, generate equivalent geometric parameters for the projection angle of each projection image. The three-dimensional coordinates of the pixels in the projected image corresponding to the detector space are calculated based on the equivalent geometric parameters of the projection angle. The normal vector of the detector plane and the normal distance from the ray source to the detector plane are calculated based on the equivalent geometric parameters of the projection angle. The corresponding projected image is weighted by cone-beam averaging based on the pixel's three-dimensional coordinates and the normal distance from the X-ray source to the detector plane to obtain a weighted projected image. The weighted projection image is filtered along the X direction of the detector to obtain a filtered weighted projection image; Based on the normal vector under the projection angle and the normal distance from the X-ray source to the detector plane, the filtered weighted projection image is accumulated by voxel projection and back projection to obtain the three-dimensional reconstructed volume data corresponding to the candidate geometric parameters.

[0023] Furthermore, the optimization module multiplies the three-dimensional position coordinates of the ray source, the three-dimensional position coordinates of the detector center, the three-dimensional displacement vector of the unit pixel in the X direction of the detector, and the three-dimensional displacement vector of the unit pixel in the Y direction of the detector in the current candidate geometric parameters with the rotation matrix corresponding to the projection angle, respectively, to obtain the equivalent geometric parameters of the projection angle. The equivalent geometric parameters of the projection angle include the three-dimensional position coordinates of the X-ray source under the projection angle, the three-dimensional position coordinates of the detector center under the projection angle, the three-dimensional displacement vector per unit pixel in the X direction of the detector under the projection angle, and the three-dimensional displacement vector per unit pixel in the Y direction of the detector under the projection angle.

[0024] Further, the optimization module calculates the three-dimensional coordinates of pixels in the projected image corresponding to the detector space based on the equivalent geometric parameters of the projection angle, including: Calculate the center pixel coordinates of the detector based on its width and height; Calculate the offset of each pixel in the projected image relative to the coordinates of the detector's center pixel; Based on the offset, the three-dimensional displacement vector of a unit pixel in the X direction of the detector under the projection angle, and the three-dimensional displacement vector of a unit pixel in the Y direction of the detector under the projection angle, the displacement vector of the pixel in the detection space is calculated. Based on the three-dimensional position coordinates of the detector center under the projection angle and the displacement vector corresponding to the pixel in the detection space, the three-dimensional coordinates of each pixel in the projected image in the detector space are obtained.

[0025] Further, the optimization module calculates the normal vector under the projection angle and the normal distance from the ray source to the detector plane based on the equivalent geometric parameters of the projection angle, including: The normal vector under the projection angle is calculated based on the three-dimensional displacement vector of a unit pixel in the X direction of the detector under the projection angle and the three-dimensional displacement vector of a unit pixel in the Y direction of the detector under the projection angle. Based on the normal vector, the three-dimensional position coordinates of the X-ray source under the projection angle, and the three-dimensional position coordinates of the detector center under the projection angle, the normal distance from the X-ray source to the detector plane is calculated.

[0026] Furthermore, the optimization module performs cone-beam weighting on the corresponding projected image based on the pixel's three-dimensional coordinates and the normal distance from the ray source to the detector plane, including: The cone-beam weighting coefficients of the corresponding pixels are calculated based on the normal distance from the X-ray source to the detector plane and the three-dimensional coordinates of the pixels. The cone-beam weighting coefficients are applied to the corresponding pixels in the projected image to obtain a weighted projected image.

[0027] Furthermore, the optimization module performs voxel projection and back projection accumulation on the filtered weighted projection image based on the normal vector at the projection angle and the normal distance from the ray source to the detector plane, including: Based on the normal vector, the three-dimensional position coordinates of the X-ray source under the projection angle, and the pre-determined center coordinates of the voxel, calculate the distance of the voxel relative to the X-ray source in the detector normal direction. If the distance of the voxel relative to the X-ray source in the detector normal direction is greater than 0, then the intersection point of the detector plane corresponding to the voxel is calculated based on the normal distance from the X-ray source to the detector plane, the distance of the voxel relative to the X-ray source in the detector normal direction, the three-dimensional position coordinates of the X-ray source under the projection angle, and the center coordinates of the voxel. The intersection points of the detector planes corresponding to the voxels are converted into continuous pixel coordinates on the filtered weighted projection image; Bilinear interpolation is performed at the continuous pixel coordinates on the filtered weighted projection image to obtain the back-projection image; The back-projected images under all projection angles are weighted and accumulated to obtain the 3D reconstructed volume data corresponding to the candidate geometric parameters.

[0028] Furthermore, the optimization module extracts sphere data from the 3D reconstructed volume data and solves for the objective function value of the error based on the sphere data, including: Based on the three-dimensional reconstructed data, threshold segmentation and three-dimensional connected component analysis are performed to calculate the reconstructed sphere center coordinates and reconstructed sphere radius of the sphere in the single-layer sphere array calibration plate. The sphericity error is calculated based on the reconstructed sphere center coordinates and reconstructed sphere radius, and the radius error is calculated based on the reconstructed sphere radius. Based on the reconstructed sphere center coordinates, establish the sphere array coordinate direction matrix and calculate the sphere center array coordinates; Calculate the standard sphere center coordinates based on the standard sphere center coordinates in the physical coordinate system of the single-layer sphere array calibration plate; The array coordinate error is calculated based on the array coordinates of the sphere center and the standard sphere center-centered coordinates, and the sphere spacing error is calculated based on the standard sphere center-centered coordinates and the reconstructed sphere center coordinates. Define the ball count penalty and parameter regularization terms; The objective function value is obtained by weighted summing of the sphericity error, radius error, array coordinate error, sphere spacing error, sphere quantity penalty term, and parameter regularization term.

[0029] Furthermore, based on the three-dimensional reconstructed data, threshold segmentation and three-dimensional connected component analysis are performed to calculate the coordinates of the reconstructed sphere center and the radius of the reconstructed sphere in the single-layer sphere array calibration plate, including: The three-dimensional reconstructed volume data is segmented according to the set segmentation threshold to obtain binary volume data representing the spherical region; Based on the binary data, a three-dimensional connected component analysis was performed to obtain multiple spherical regions; Based on the center coordinates of each voxel extracted from the sphere region and the corresponding 3D reconstructed volume data, the coordinates of the reconstructed sphere center are calculated by gray-scale weighting. The radius of the reconstructed sphere is obtained by performing least-squares fitting on the boundary points of the sphere region.

[0030] Furthermore, the optimization module establishes a sphere array coordinate direction matrix based on the reconstructed sphere center coordinates, including: Obtain the coordinates of all reconstructed sphere centers to form a set of sphere centers; Calculate the sphere center covariance matrix based on the sphere center set, perform eigenvalue decomposition on the sphere center covariance matrix, and take the unit eigenvectors corresponding to the two largest eigenvalues ​​as the row and column directions of the sphere array coordinate system; The height direction is obtained by cross-multiplying the unit eigenvectors of the row and column directions. The row, column, and height directions together form the coordinate direction matrix of the spherical array.

[0031] Furthermore, the sphericity error is the sum of the squares of the normalized spherical residuals of the sphere, and the normalized spherical residual of a single sphere is the ratio of the root mean square residual of the distance from the boundary point of the sphere region to the center of the reconstructed sphere relative to the radius of the reconstructed sphere to the radius of the reconstructed sphere. The radius error is the mean square error between the reconstructed sphere radius and the standard sphere radius of all spheres; The sphere center array coordinates are the product of the transpose of the sphere array coordinate direction matrix and the reconstructed sphere center coordinates and the difference between the reconstructed sphere center coordinates and the mean of the reconstructed sphere center coordinates; The standard sphere center-centered coordinates are the difference between the standard sphere center coordinates and the mean of the standard sphere center coordinates in the physical coordinate system of the single-layer sphere array calibration plate. The array coordinate error is the mean square error of the difference between the array coordinates of the center of all spheres and the standard sphere center-centered coordinates; the sphere spacing error is the mean square error of the difference between the reconstructed sphere center spacing of all two adjacent spheres and the standard sphere center spacing. The ball count penalty term is the absolute value of the difference between the number of balls detected and the actual number of balls under the current candidate geometric parameters; the parameter regularization term is the sum of squared deviations of the changes between the initial geometric parameters and the current candidate parameters.

[0032] The method and apparatus for calibrating geometric parameters of a CT system based on reconstruction feedback optimization provided by this invention have at least the following beneficial effects: (1) The optimization of geometric parameters is directly based on the three-dimensional reconstruction results, so that the calibration target no longer depends solely on the two-dimensional projection error, but rather the geometric parameters are evaluated jointly by the center position of the sphere array, the distance between the spheres, the sphericity and the radius consistency, thereby improving the accuracy, robustness and engineering applicability of CT geometric calibration; (2) The geometric parameters to be calibrated are clearly defined as the position of the X-ray source, the center position of the detector, the X-axis vector of the detector, and the Y-axis vector of the detector, so as to avoid treating the rotation angle, the geometry of each angle, the normal vector, the voxel projection point and the PCA principal direction as independent optimization parameters, thereby reducing the parameter dimension and parameter coupling. (3) Based on the CT system geometric coordinate system established on the rotation axis, it can clearly describe the spatial orientation of the X-ray source, detector center and detector pixel direction in the same coordinate system, so that the geometric model directly corresponds to the X-ray calculation process in CT reconstruction. (4) A single-layer M×N sphere array calibration plate is used. Geometric parameters are evaluated by the sphere center, sphere radius, sphericity and array spatial distribution in the reconstructed volume, so that the calibration target is directly oriented towards the final three-dimensional reconstruction quality, rather than relying solely on the two-dimensional projection domain error. (5) First, use PCA to construct the sphere array coordinate direction matrix from the reconstructed sphere center, then transform the reconstructed sphere center to the array coordinate system, and subtract it from the coordinates of the centered standard sphere array. This can reduce the influence of the calibration plate placement posture on the error evaluation. (6) The array coordinate error, sphere spacing error, sphericity error, radius error, sphere number penalty and parameter regularization term are combined into a unified objective function, so that the geometric parameter optimization takes into account the overall positional relationship, local sphere shape and equipment physical constraints, thus improving the stability of the calibration results. Attached Figure Description

[0033] Figure 1 This is a flowchart of one embodiment of the geometric parameter calibration method for a CT system based on reconstruction feedback optimization provided by the present invention.

[0034] Figure 2 This is a schematic diagram of an embodiment of the CT system geometric coordinate system in the CT system geometric parameter calibration method based on reconstruction feedback optimization provided by the present invention.

[0035] Figure 3 This is a schematic diagram of one embodiment of the geometric parameter calibration device for a CT system based on reconstruction feedback optimization provided by the present invention. Detailed Implementation

[0036] To better understand the above technical solutions, the following will provide a detailed explanation of the technical solutions in conjunction with the accompanying drawings and specific implementation methods.

[0037] refer to Figure 1 In some embodiments, a method for calibrating geometric parameters of a CT system based on reconstruction feedback optimization is provided, including: S1. Establish the geometric coordinate system of the CT system and determine the geometric parameters to be calibrated; S2. Set up a single-layer ball array calibration plate; S3. Control the CT system to acquire projection images of the single-layer spherical array calibration plate at different projection angles according to the scanning trajectory, and form a set of projection images; S4. Iteratively execute the following operations until the stopping condition is met: Take the position of each particle in the particle swarm as a set of candidate geometric parameters. Under the constraints of the geometric coordinate system of the CT system, perform three-dimensional reconstruction based on the set of projected images and the corresponding candidate geometric parameters to obtain the three-dimensional reconstructed volume data corresponding to each set of candidate geometric parameters; extract sphere data based on the three-dimensional reconstructed volume data; solve the objective function value with respect to the error based on the sphere data; and update the position of the particles. S5. The globally optimal solution obtained when the stopping condition is met is used as the final calibration geometric parameters.

[0038] Specifically, in step S1, refer to Figure 2 In the established CT system geometric coordinate system, the origin O is the intersection of the plane containing the rotation axis and the center of rotation, the Z-axis is the rotation axis, and the X-axis and Y-axis are located in a plane perpendicular to the rotation axis and form a right-handed coordinate system with the Z-axis. Since the CT system geometric coordinate system is established on the rotation axis, when the rotating platform rotates, the measured object or the single-layer sphere array calibration plate changes angle relative to this coordinate system, and this angle is the projection angle. In the equivalent reconstruction description, the measured object is considered to be stationary, while the X-ray source and detector rotate relative to each other around the Z-axis with the projection angle.

[0039] In the geometric coordinate system of the CT system, the geometric parameter Θ to be calibrated d This includes the three-dimensional position coordinates S of the X-ray source at the initial projection angle, the three-dimensional position coordinates D of the detector center at the initial projection angle, the three-dimensional displacement vector U per unit pixel in the X direction of the detector at the initial projection angle, and the three-dimensional displacement vector V per unit pixel in the Y direction of the detector at the initial projection angle. (1) U and V describe the detector's attitude, pixel orientation, and pixel size. Therefore, the detector's translation, in-plane rotation, pitch, roll, and pixel orientation errors can all be represented by D, U, and V. Expanding these four three-dimensional vectors, the geometric parameter vector to be calibrated can be written as: (2) Wherein, Sx, Sy, and Sz represent the x, y, and z coordinates of the X-ray source under the initial projection angle; Dx, Dy, and Dz represent the x, y, and z coordinates of the detector center under the initial projection angle; Ux, Uy, and Uz represent the x, y, and z coordinates of the three-dimensional displacement vector per unit pixel in the X direction of the detector under the initial projection angle; and Vx, Vy, and Vz represent the x, y, and z coordinates of the three-dimensional displacement vector per unit pixel in the Y direction of the detector under the initial projection angle.

[0040] Furthermore, in step S2, the single-layer ball array calibration plate includes a substrate and balls arranged in an array on the substrate.

[0041] Specifically, to ensure that geometric errors can form stable and quantifiable feedback in the reconstructed volume, this embodiment employs a single-layer sphere array calibration plate. This single-layer sphere array calibration plate comprises a low-absorptivity substrate and multiple high-absorptivity spheres. The low-absorptivity substrate is made of engineering plastic, while the high-absorptivity spheres are made of steel. Due to the significant difference in X-ray attenuation coefficients between the spheres and the substrate, high contrast is achieved in both the projected image and the reconstructed volume data.

[0042] Suppose that a single-layer sphere array calibration plate has M rows and N columns of spheres, with a total of K spheres: K = M × N; (3) A physical coordinate system is established using a single-layer sphere array calibration plate. The origin of the physical coordinate system is taken as the center of a certain reference sphere, the column direction is the X-axis, the row direction is the Y-axis, and the direction perpendicular to the plane of the single-layer sphere array calibration plate is the Z-axis. Since the single-layer sphere array calibration plate is a single-layer sphere array, the Z-coordinate of all standard sphere centers is set to 0.

[0043] Let the column spacing of the single-layer ball array calibration plate be a. x The line spacing is a y The standard coordinates of the center of the sphere in the m-th row and n-th column in the physical coordinate system of the calibration plate are denoted as Q. m,n Where m = 0, 1, ..., M-1, n = 0, 1, ..., N-1, then: (4) To avoid directional ambiguity when constructing the coordinate orientation matrix of the spherical array using subsequent PCA (Principal Component Analysis), M and N are chosen to be unequal. This ensures that the row and column directions of the array can be stably determined based on the distribution of the sphere centers after reconstruction.

[0044] In step S3, the CT system is controlled to acquire multiple projection images of a single-layer spherical array calibration plate at different projection angles according to the scanning trajectory, forming a projection image set.

[0045] Further, in step S4, a particle algorithm is used to optimize the geometric parameters. In each iteration, the position of each particle in the particle swarm is used as a set of candidate geometric parameters. Each projection image in the projection image set corresponds to a projection angle. Under the constraints of the CT system's geometric coordinate system, 3D reconstruction is performed based on the projection image set and the corresponding candidate geometric parameters to obtain the 3D reconstructed volume data corresponding to each set of candidate geometric parameters. Sphere data is extracted from the 3D reconstructed volume data, and the objective function value with respect to error is solved based on the sphere data. The particle position is updated based on the individual historical optimal solution and the historical global optimal solution. This process is repeated until a stopping condition is met. This stopping condition allows the change in the objective function value of the global optimal solution to be less than a threshold, or the preset number of iterations to be reached.

[0046] The 3D reconstruction based on the set of projected images and corresponding candidate geometric parameters includes: S41. Based on the candidate geometric parameters and the projection angle of each projection image in the projection image set, generate equivalent geometric parameters for the projection angle of each projection image. S42. Calculate the three-dimensional coordinates of the pixels in the projected image corresponding to the detector space based on the equivalent geometric parameters of the projection angle; S43. Calculate the normal vector of the detector plane under the projection angle and the normal distance from the ray source to the detector plane based on the equivalent geometric parameters of the projection angle. S44. Based on the pixel's three-dimensional coordinates and the normal distance from the ray source to the detector plane, the corresponding projection image is subjected to cone-beam weighting to obtain a weighted projection image; S45. Filter the weighted projection image along the X direction of the detector to obtain a filtered weighted projection image; S46. Based on the normal vector under the projection angle and the normal distance from the X-ray source to the detector plane, perform voxel projection and back projection accumulation on the filtered weighted projection image to obtain the three-dimensional reconstructed volume data corresponding to the candidate geometric parameters.

[0047] Further, in step S41, the three-dimensional position coordinates of the ray source, the three-dimensional position coordinates of the detector center, the three-dimensional displacement vector of the unit pixel in the X direction of the detector, and the three-dimensional displacement vector of the unit pixel in the Y direction of the detector in the candidate geometric parameters are multiplied by the rotation matrix corresponding to the projection angle to obtain the equivalent geometric parameters of the projection angle. The equivalent geometric parameters of the projection angle include the three-dimensional position coordinates of the X-ray source under the projection angle, the three-dimensional position coordinates of the detector center under the projection angle, the three-dimensional displacement vector per unit pixel in the X direction of the detector under the projection angle, and the three-dimensional displacement vector per unit pixel in the Y direction of the detector under the projection angle.

[0048] Let θ be the projection angle corresponding to the q-th projected image. qWhere q = 1, 2, 3...Q, and Q is the total number of projected images. Let R be... q The projection angle θ is the rotation about the Z-axis of the world coordinate system. q The rotation matrix is ​​used to generate the equivalent geometry at the projection angle, based on the definition of the scanning direction, as shown below: S q =R q S; D q =R q D; U q =R q U; V q =R q V; (5) Among them, S q Represents the projection angle θ q The three-dimensional position coordinates of the ray source are shown below, D. q U represents the three-dimensional position coordinates of the detector center under the projection angle. q Let θ represent the 3D displacement vector per unit pixel in the X direction of the detector under the projection angle, and Vq represent the 3D displacement vector per unit pixel in the Y direction of the detector under the projection angle. Therefore, the projection angle θ... q The equivalent geometric parameter Gq of the projection angle under the given conditions can be expressed as: (6) In step S42, calculating the three-dimensional coordinates of pixels in the projected image corresponding to the detector space based on the equivalent geometric parameters of the projection angle includes: S421. Calculate the coordinates of the center pixel of the detector based on the width and height of the detector; S422. Calculate the offset of each pixel in the projected image relative to the coordinates of the detector center pixel; S423. Based on the offset, the three-dimensional displacement vector of a unit pixel in the X direction of the detector under the projection angle, and the three-dimensional displacement vector of a unit pixel in the Y direction of the detector under the projection angle, calculate the displacement vector of the pixel in the detection space. S424 obtains the three-dimensional coordinates of each pixel in the projected image in the detector space based on the three-dimensional position coordinates of the detector center under the projection angle and the displacement vector corresponding to the pixel in the detection space.

[0049] Specifically, let the detector image width be N. u The height is N v The detector pixel coordinates are (i, j), where i = 0, 1, 2, ..., N. u -1, j=0,1,2……N v-1, the pixel coordinates (i0, j0) corresponding to the center of the detector are denoted as: (7) Then the projection angle θ q The three-dimensional coordinates Pq(i,j) of the pixel in the detector space corresponding to the pixel in the projected image are: (8) Where (i-i0) and (j-j0) represent the offset of each pixel in the projected image relative to the coordinates of the detector center pixel, and (i-i0)Uq and (j-j0)Vq represent the displacement vector of the pixel in the detection space.

[0050] Further, in step S43, calculating the normal vector under the projection angle and the normal distance from the ray source to the detector plane based on the equivalent geometric parameters of the projection angle includes: S431. Based on the three-dimensional displacement vector U of the unit pixel in the X direction of the detector under the projection angle. q And the three-dimensional displacement vector V per unit pixel in the Y direction of the detector under the projection angle. q Calculate the normal vector under the projection angle; S432, Based on the normal vector and projection angle, the three-dimensional position coordinates of the ray source S q and the three-dimensional position coordinates D of the detector center under the projection angle. q Calculate the normal distance from the X-ray source to the detector plane.

[0051] Specifically, for the current candidate geometric parameters, the equivalent geometric parameters Gq of the projection angle are generated, and the normal vector of the detector plane under the projection angle is defined as n. q To ensure that the normal vector direction points from the ray source to the detector, if the calculated normal vector direction is perpendicular to D... q -S q Conversely, take the opposite direction. This can be written as: (9) And make it satisfy: (10) The normal distance Lq from the X-ray source to the detector plane is: (11) Further, in step S44, cone-beam weighting is performed on the corresponding projected image based on the pixel's three-dimensional coordinates and the normal distance from the ray source to the detector plane, including: S441. Calculate the cone-beam weighting coefficient of the corresponding pixel based on the normal distance from the X-ray source to the detector plane and the three-dimensional coordinates of the pixel. S442. Apply the cone-beam weighting coefficients to the corresponding pixels in the projected image to obtain a weighted projected image.

[0052] Specifically, for each pixel in the projected image, the cone-beam weighting coefficient of that pixel is calculated from its corresponding three-dimensional pixel coordinates Pq(i,j) in the detector space: (12) Among them, w q (i,j) are the cone-beam weighting coefficients; Applying cone-beam weighting to the entire projected image yields a weighted projected image: (13) in, p represents a weighted projection image. q (i,j) represents the original projected image. Since the original projected image is weighted, the resulting weighted projected image is still a two-dimensional image.

[0053] Further, in step S45, for each projection angle, the weighted projection image is... One-dimensional ramp filtering is performed row by row along the X-axis of the detector, i.e., along the pixel column direction. Let the filter kernel be h, then the filtered weighted projection image... Its spatial domain is expressed as: (14) Where m is the column index of the column participating in the convolution within the same row. This filtering process can also be implemented in the frequency domain. Let... This represents a one-dimensional Fourier transform along the i-direction. If we represent a frequency variable, then: (15) in, This is the frequency response of the Ramp filter. This indicates taking the filtered weighted projection image. All columns in Indicates taking the weighted projection image All columns.

[0054] Further, in step S46, based on the normal vector at the projection angle and the normal distance from the ray source to the detector plane, voxel projection and back projection accumulation are performed on the filtered weighted projection image, including: S461. Based on the normal vector and the three-dimensional position coordinates of the X-ray source under the projection angle, as well as the predetermined center coordinates of the voxel, calculate the distance of the voxel relative to the X-ray source in the detector normal direction. S462. If the distance of the voxel relative to the X-ray source in the detector normal direction is greater than 0, then calculate the intersection point of the detector plane corresponding to the voxel based on the normal distance from the X-ray source to the detector plane, the distance of the voxel relative to the X-ray source in the detector normal direction, the three-dimensional position coordinates of the X-ray source under the projection angle, and the center coordinates of the voxel. S463. Convert the intersection points of the detector planes corresponding to the voxels into continuous pixel coordinates on the filtered weighted projection image; S464. Perform bilinear interpolation at the continuous pixel coordinates on the filtered weighted projection image to obtain the back-projection image; S465. Weighted summation of the back-projected images under all projection angles to obtain the 3D reconstructed volume data corresponding to the candidate geometric parameters.

[0055] Specifically, the center coordinates of a voxel are pre-determined sampling coordinates based on the preset reconstruction space range, voxel size, and three-dimensional voxel mesh index. Let the center coordinates of the voxel be X, which is represented in the geometric coordinate system of the CT system. For the q-th projection image, the ray originates from the source point S. q After passing through voxel X, it intersects with the detector plane; the intersection point is denoted as H. q (x).

[0056] Define process variables: (16) Where, d q (X) represents the distance of the voxel relative to the X-ray source in the detector normal direction. If d q If (X)≤0, then the voxel is located in the invalid region of the current projective geometry and can be excluded from the back projection at that angle.

[0057] Otherwise, the intersection point H of the detector plane corresponding to the voxel q (X) is: (17) Due to the intersection point H of the detector plane corresponding to the voxel q (X) is a continuous spatial point, which needs to be converted into continuous pixel coordinates in a filtered weighted projection image. Due to the intersection relationship, we have: (18) Where, α q (X) and β q (X) represent the continuous column coordinates and continuous row coordinates of the voxel projected onto the q-th filtered weighted projection image, respectively.

[0058] remember: (19) And define: (20) Among them, B q U is a three-dimensional displacement vector per unit pixel. q and V q The matrix formed by the inner product, b q (X) is H q (X)-D q respectively with U q and V q The vector formed by the inner product of .

[0059] The offset ΔH of the intersection point in the detector pixel coordinate system is: ;(twenty one)

[0060] Subsequently, in the filtered weighted projection image Upper position (α) q (X) and β q (X) is subjected to bilinear interpolation to obtain the back-projected image. : ;(twenty two) Let r q This represents the distance from the source point to the rotation center O along the detector's normal direction: ;(twenty three) Then the three-dimensional reconstructed data It can be obtained by summing up all the projection angles: ;(twenty four) in, This represents the angular interval between adjacent projection angles. From the above formula, it can be seen that, given a set of projected images, the 3D reconstructed volume data... Completely determined by geometric parameters The decision is made based on the combination of geometric parameters. Therefore, different combinations of geometric parameters will produce different reconstructed sphere array shapes and spatial positions.

[0061] The 3D reconstructed volume data includes the reconstruction scalar value at the center coordinate X of each voxel in the pre-defined 3D voxel grid. This reconstruction scalar value can characterize the attenuation coefficient or the reconstruction grayscale.

[0062] Further, in step S4, sphere data is extracted from the 3D reconstructed volume data, and the objective function value with respect to the error is solved based on the sphere data, including: S47. Based on the three-dimensional reconstructed data, perform threshold segmentation and three-dimensional connected component analysis, calculate the reconstructed sphere center coordinates and reconstructed sphere radius of the sphere in the single-layer sphere array calibration plate, and calculate the sphericity error based on the reconstructed sphere center coordinates and reconstructed sphere radius, and calculate the radius error based on the reconstructed sphere radius; S48. Based on the reconstructed sphere center coordinates, establish the sphere array coordinate direction matrix and calculate the sphere center array coordinates; S49. Calculate the standard sphere center coordinates based on the standard sphere center coordinates in the physical coordinate system of the single-layer sphere array calibration plate. S410. Calculate the array coordinate error based on the array coordinates of the sphere center and the standard sphere center-centered coordinates, and calculate the sphere spacing error based on the standard sphere center-centered coordinates and the reconstructed sphere center coordinates. S411, Define the ball quantity penalty term and parameter regularization term; S412. The objective function value is obtained by weighted summation of the sphericity error, radius error, array coordinate error, sphere spacing error, sphere quantity penalty term, and parameter regularization term.

[0063] Specifically, in step S47, the sphere data includes the reconstructed sphere center coordinates and the reconstructed sphere radius. Based on the three-dimensional reconstructed volume data, threshold segmentation and three-dimensional connected component analysis are performed to calculate the reconstructed sphere center coordinates and the reconstructed sphere radius in the single-layer sphere array calibration plate, including: S471. The three-dimensional reconstructed volume data is segmented according to the set segmentation threshold to obtain binary volume data representing the spherical region; S472. Perform three-dimensional connected component analysis based on the binary data to obtain multiple spherical regions; S473. Based on the center coordinates of each voxel in the sphere region and the corresponding 3D reconstructed volume data, calculate the coordinates of the reconstructed sphere center using gray-scale weighting. S474. Perform least-squares fitting on the boundary points of the sphere region to obtain the radius of the reconstructed sphere.

[0064] Specifically, after obtaining the 3D reconstructed volume data corresponding to the current candidate geometric parameters, 3D feature extraction is performed on the single-layer sphere array calibration plate. First, threshold segmentation is performed on the 3D reconstructed volume data to obtain the high-absorption sphere region. Let the segmentation threshold be T, and the binary volume data be... ,but: (25) Where X represents the center coordinates of any voxel in the 3D reconstructed volume. The threshold T is determined based on the peak value of the gray-level histogram, and the valley between the gray-level peak of the background region and the gray-level peak of the sphere is selected as the threshold.

[0065] Three-dimensional connected component analysis is performed on binary data. Specifically, the voxels corresponding to the binary data that meet the threshold conditions are marked as foreground voxels. Connected component marking is performed on adjacent foreground voxels to obtain several connected regions. Candidate calibration spherical regions are selected based on the volume of the connected regions, the size of the bounding box, the gray mean, and the shape features. Noise points, adhesion artifacts, and non-spherical regions are removed to obtain multiple spherical regions, each of which includes a sphere.

[0066] For the k-th sphere region Calculate its gray-weighted sphere center. Let X be... l Indicates the ball area Inner l Individual center coordinates, l =1,2,……N k N k Let Ck be the number of voxels within the sphere region. Then, the reconstructed center coordinates Ck of the k-th sphere can be expressed as: (26) in, This represents the 3D reconstructed volume data corresponding to the voxel.

[0067] To further evaluate the sphere's shape, we can extract the set of boundary points for the k-th sphere region. Let the l-th boundary point in the boundary point set be Y. k,1 The number of boundary points is M k The radius r of the reconstructed sphere is obtained by performing least-squares sphere fitting on the boundary points. k : (27) Furthermore, after obtaining the coordinates of the reconstructed sphere center and the reconstructed sphere radius, the sphericity error and radius error are calculated. The sphericity error is the mean of the sum of squares of the normalized spherical residuals of the sphere. The normalized spherical residual of a single sphere is the ratio of the root mean square residual of the distance from the boundary point of the sphere region to the reconstructed sphere center relative to the reconstructed sphere radius to the reconstructed sphere radius.

[0068] Normalized spherical residual of the k-th sphere for: (28) When the geometric parameters are accurate, the reconstructed sphere approximates an ideal sphere, and the normalized spherical residuals are accurate. The reconstructed sphere may exhibit stretching, ghosting, or ellipsoidization when the geometric parameters are inaccurate, resulting in a smaller normalized spherical residual. Increase.

[0069] Sphericity error of candidate geometric parameter Θ for: (29) Where K is the number of spheres in the single-layer sphere array calibration plate.

[0070] radius error The mean square error between the reconstructed sphere radius and the standard sphere radius for all spheres: (30) Where r0 is the standard sphere radius, that is, the actual radius of the sphere. If different calibrated spheres have different radii, r0 can be replaced with the standard radius corresponding to each sphere.

[0071] Further, in step S48, establishing the sphere array coordinate direction matrix based on the reconstructed sphere center coordinates includes: S481. Obtain the coordinates of all reconstructed sphere centers to form a set of sphere centers; S482. Calculate the sphere center covariance matrix based on the sphere center set, perform eigenvalue decomposition on the sphere center covariance matrix, and take the unit eigenvectors corresponding to the two largest eigenvalues ​​as the row direction and column direction of the sphere array coordinate direction matrix; S483. Cross-multiply the unit eigenvectors of the row direction and column direction to obtain the height direction. The row direction, column direction and height direction form the coordinate direction matrix of the spherical array.

[0072] Specifically, since the initial orientation of the single-layer spherical array calibration plate placed on the rotating platform may involve translation and rotation, the reconstructed sphere center coordinates C extracted from the CT system's geometric coordinate system cannot be directly obtained. k With standard sphere center coordinates Q m,n To eliminate the influence of the orientation of the single-layer sphere array calibration board, this embodiment uses the extracted set of effective sphere centers to construct the coordinate direction matrix of the reconstructed sphere array.

[0073] The set of centers of the sphere is denoted as: (31) in, Let K be the set of sphere centers, and K be the actual number of spheres in the single-layer sphere array calibration plate. If the number of detected spheres is not equal to K, then the set of geometric parameters corresponds to a larger penalty.

[0074] Calculate the mean value of the ball's center : (32) Calculate the covariance matrix of the sphere center : (33) Covariance matrix of the center of the sphere Eigenvalue decomposition yields two main directions. Since the single-layer sphere array calibration plate is a single-layer M×N sphere array, with the sphere centers primarily distributed within a single plane, the unit eigenvectors corresponding to the two largest eigenvalues ​​are taken as the column and row directions in the sphere array coordinate direction matrix, denoted as ei, respectively. x and e y And order: (34) Obtain height direction e z.

[0075] Therefore, the coordinate direction matrix R of the spherical array is constructed. rec : (35) Wherein, the spherical array coordinate direction matrix R rec The three columns represent the x, y, and z directions of the reconstructed spherical array coordinate system. If there is ambiguity in the orientation of the spherical array, it can be corrected based on the differences in the number of M and N, the differences in row and column spacing, or the position of the reference sphere. x and e y The direction.

[0076] After obtaining the coordinate direction matrix R of the spherical array rec Then, each reconstructed sphere center is transformed from the CT system geometric coordinate system to the reconstructed sphere array coordinate system. The sphere center array coordinates of the k-th sphere center are the product of the sphere array coordinate direction matrix and the difference between the reconstructed sphere center coordinates and the mean of the reconstructed sphere center coordinates.

[0077] The coordinates Z of the center array of the k-th sphere are... k The calculation method is as follows: (36) in, To reconstruct the mean coordinates of the sphere's center.

[0078] In step S49, the standard sphere center coordinates in the physical coordinate system of the single-layer sphere array calibration plate are renumbered according to the sphere array coordinate direction, and the standard sphere center-centered coordinates are calculated. The standard sphere center-centered coordinates are the difference between the renumbered standard sphere center coordinates and the mean of the standard sphere center coordinates.

[0079] Specifically, for the standard sphere center coordinates Q in the m-th row and n-th column of the physical coordinate system of the single-layer sphere array calibration plate... m,n Where m = 0, 1, ..., M−1, n = 0, 1, ..., N-1, according to row precedence rules, let k = mN+n+1, and denote Q. k =Q m,n After completing the array orientation correction, the array coordinates of the reconstructed sphere centers are clustered and sorted according to their row and column positions, and k is numbered according to the same rules to ensure that the reconstructed sphere centers correspond one-to-one with the standard sphere centers.

[0080] The standard sphere center coordinates represented by the above one-dimensional numbering are: Its mean for: (37) Standard sphere center-centered coordinates for: (38) At this moment, the coordinates Z of the sphere center array k and standard sphere center-centered coordinates Each sphere lies within its own array coordinate system, which can be used to construct the objective function. In other words, this embodiment does not directly compare the centers of the spheres in the CT system's geometric coordinate system, but first uses PCA to obtain the sphere array coordinate direction matrix R. rec Then use the spherical array coordinate direction matrix R rec The reconstructed sphere center is converted to array coordinates, and finally the difference between the array coordinates and the standard sphere center is calculated.

[0081] Further, in step S410, the array coordinate error The mean square error of the difference between the array coordinates of the centers of all spheres and the standard sphere-centered coordinates: (39) The ball spacing error The mean square error of the difference between the reconstructed center-to-center distance and the standard center-to-center distance for all adjacent spheres: (40) Where N is the selected set of ball pairs, which can be adjacent ball pairs, adjacent ball pairs in the same row, or adjacent ball pairs in the same column. Standard ball center distance. for That is, the actual distance between the centers of two adjacent spheres in a single-layer sphere array calibration plate, reconstructing the sphere center spacing. for .

[0082] In step S411, the ball quantity penalty item This represents the absolute value of the difference between the number of detected spheres and the actual number of spheres under the current candidate geometric parameters. (41) in, K represents the number of spheres detected under the current candidate geometric parameters, and K represents the actual number of spheres in the single-layer sphere array calibration plate.

[0083] Parameter regularization terms The sum of squared deviations between the initial geometric parameters and the current candidate parameters: (42) in, These are the initial geometric parameters. Let θj and θ be the allowable scale of variation for the j-th parameter. j,0 These are the current candidate geometric parameters. and initial geometric parameters The j-th parameter.

[0084] In summary, in step S412, the objective function for: (43) Where λ1 to λ6 are the weighting coefficients of each error term.

[0085] Among them, the sphere spacing error is used to evaluate whether the spatial relative relationship between the centers of the reconstructed spheres is consistent with the standard array; the sphericity error is used to evaluate whether the reconstructed sphere shape maintains the ideal spherical shape; the radius error is used to evaluate whether the size of the reconstructed sphere is consistent with the standard radius; the sphere number penalty term is used to constrain the actual number of detected spheres to be consistent with the standard number of spheres; and the parameter regularization term is used to limit the optimization parameters from deviating too far from the initial geometry.

[0086] The smaller the objective function, the closer the reconstructed spherical array is to a standard spherical array under the current geometric parameters, the closer the spherical shape is to a real sphere, and the closer the system's geometric model is to the actual CT imaging state. Therefore, the objective of step S4 is to find the geometric parameters that minimize the objective function.

[0087] After calculating the objective function value, let the best historical position of the i-th particle be... The historical global optimal position is Then the particle swarm update formula is: (44) (45) in, Here, c1 is the inertia weight, c2 is the individual learning factor, and r1 and r2 are random numbers between 0 and 1. Indicates the particle refresh rate. Indicates the current updated position of the particle. This represents the particle's previous position. To balance global search and local convergence, the inertia weight can be gradually decreased with the number of iterations: (46) Among them, T max The maximum number of iterations, and These are the initial inertia weight and the final inertia weight, respectively.

[0088] When a particle's position exceeds the allowed range of parameters, it is either truncated at the boundary or subjected to velocity bounce.

[0089] In step S5, when the change in the global optimal objective function is less than a threshold in several consecutive iterations, or when the maximum number of iterations is reached, the optimization stops, and the global optimal solution is output. (47) Finally, using the final calibrated geometric parameters The corrected geometric model is generated and used for subsequent CT reconstruction of the workpiece under test.

[0090] refer to Figure 3 In some embodiments, a CT system geometric parameter calibration device based on reconstruction feedback optimization applied to the above method is provided, comprising: The coordinate establishment module 201 is used to establish the geometric coordinate system of the CT system and determine the geometric parameters to be calibrated. Image acquisition module 202 is used to control the CT system to acquire projection images of a single-layer spherical array calibration plate at different projection angles according to the scanning trajectory, and form a set of projection images; The optimization module 203 is used to iteratively perform the following operations until the stopping condition is met: taking the position of each particle in the particle swarm as a set of candidate geometric parameters, and under the constraints of the geometric coordinate system of the CT system, performing three-dimensional reconstruction based on the set of projected images and the corresponding candidate geometric parameters to obtain the three-dimensional reconstructed volume data corresponding to each set of candidate geometric parameters; extracting sphere data based on the three-dimensional reconstructed volume data, solving the objective function value with respect to the error based on the sphere data, and updating the position of the particles; The parameter determination module 204 is used to take the globally optimal solution obtained when the stopping condition is met as the final calibration geometric parameters.

[0091] Furthermore, in the geometric coordinate system of the CT system, the origin is the intersection of the plane containing the rotation axis and the center of rotation, the Z-axis is the rotation axis, and the X-axis and Y-axis are located in a plane perpendicular to the rotation axis and form a right-handed coordinate system with the Z-axis; In the geometric coordinate system of the CT system, the geometric parameters to be calibrated include the three-dimensional position coordinates of the X-ray source under the initial projection angle, the three-dimensional position coordinates of the detector center under the initial projection angle, the three-dimensional displacement vector per unit pixel in the X direction of the detector under the initial projection angle, and the three-dimensional displacement vector per unit pixel in the Y direction of the detector under the initial projection angle.

[0092] Furthermore, the single-layer sphere array calibration plate includes a substrate and spheres arranged in an array on the substrate.

[0093] Furthermore, the optimization module 203 performs 3D reconstruction based on the set of projected images and the corresponding candidate geometric parameters, including: Based on the candidate geometric parameters and the projection angle of each projection image in the projection image set, generate equivalent geometric parameters for the projection angle of each projection image. The three-dimensional coordinates of the pixels in the projected image corresponding to the detector space are calculated based on the equivalent geometric parameters of the projection angle. The normal vector of the detector plane and the normal distance from the ray source to the detector plane are calculated based on the equivalent geometric parameters of the projection angle. The corresponding projected image is weighted by cone-beam averaging based on the pixel's three-dimensional coordinates and the normal distance from the X-ray source to the detector plane to obtain a weighted projected image. The weighted projection image is filtered along the X direction of the detector to obtain a filtered weighted projection image; Based on the normal vector under the projection angle and the normal distance from the X-ray source to the detector plane, the filtered weighted projection image is accumulated by voxel projection and back projection to obtain the three-dimensional reconstructed volume data corresponding to the candidate geometric parameters.

[0094] Furthermore, the optimization module 203 multiplies the three-dimensional position coordinates of the ray source, the three-dimensional position coordinates of the detector center, the three-dimensional displacement vector of the unit pixel in the X direction of the detector, and the three-dimensional displacement vector of the unit pixel in the Y direction of the detector in the current candidate geometric parameters with the rotation matrix corresponding to the projection angle, respectively, to obtain the equivalent geometric parameters of the projection angle. The equivalent geometric parameters of the projection angle include the three-dimensional position coordinates of the X-ray source under the projection angle, the three-dimensional position coordinates of the detector center under the projection angle, the three-dimensional displacement vector per unit pixel in the X direction of the detector under the projection angle, and the three-dimensional displacement vector per unit pixel in the Y direction of the detector under the projection angle.

[0095] Further, the optimization module 203 calculates the three-dimensional coordinates of pixels in the projected image corresponding to the detector space based on the equivalent geometric parameters of the projection angle, including: Calculate the center pixel coordinates of the detector based on its width and height; Calculate the offset of each pixel in the projected image relative to the coordinates of the detector's center pixel; Based on the offset, the three-dimensional displacement vector of a unit pixel in the X direction of the detector under the projection angle, and the three-dimensional displacement vector of a unit pixel in the Y direction of the detector under the projection angle, the displacement vector of the pixel in the detection space is calculated. Based on the three-dimensional position coordinates of the detector center under the projection angle and the displacement vector corresponding to the pixel in the detection space, the three-dimensional coordinates of each pixel in the projected image in the detector space are obtained.

[0096] Further, the optimization module 203 calculates the normal vector under the projection angle and the normal distance from the ray source to the detector plane based on the equivalent geometric parameters of the projection angle, including: The normal vector under the projection angle is calculated based on the three-dimensional displacement vector of a unit pixel in the X direction of the detector under the projection angle and the three-dimensional displacement vector of a unit pixel in the Y direction of the detector under the projection angle. Based on the normal vector, the three-dimensional position coordinates of the X-ray source under the projection angle, and the three-dimensional position coordinates of the detector center under the projection angle, the normal distance from the X-ray source to the detector plane is calculated.

[0097] Furthermore, the optimization module 203 performs cone-beam weighting on the corresponding projected image based on the pixel's three-dimensional coordinates and the normal distance from the ray source to the detector plane, including: The cone-beam weighting coefficients of the corresponding pixels are calculated based on the normal distance from the X-ray source to the detector plane and the three-dimensional coordinates of the pixels. The cone-beam weighting coefficients are applied to the corresponding pixels in the projected image to obtain a weighted projected image.

[0098] Furthermore, the optimization module 203 performs voxel projection and back projection accumulation on the filtered weighted projection image based on the normal vector at the projection angle and the normal distance from the X-ray source to the detector plane, including: Based on the normal vector, the three-dimensional position coordinates of the X-ray source under the projection angle, and the pre-determined center coordinates of the voxel, calculate the distance of the voxel relative to the X-ray source in the detector normal direction. If the distance of the voxel relative to the X-ray source in the detector normal direction is greater than 0, then the intersection point of the detector plane corresponding to the voxel is calculated based on the normal distance from the X-ray source to the detector plane, the distance of the voxel relative to the X-ray source in the detector normal direction, the three-dimensional position coordinates of the X-ray source under the projection angle, and the center coordinates of the voxel. The intersection points of the detector planes corresponding to the voxels are converted into continuous pixel coordinates on the filtered weighted projection image; Bilinear interpolation is performed at the continuous pixel coordinates on the filtered weighted projection image to obtain the back-projection image; The back-projected images under all projection angles are weighted and accumulated to obtain the 3D reconstructed volume data corresponding to the candidate geometric parameters.

[0099] Furthermore, the optimization module 203 extracts sphere data from the 3D reconstructed volume data and solves for the objective function value of the error based on the sphere data, including: Based on the three-dimensional reconstructed data, threshold segmentation and three-dimensional connected component analysis are performed to calculate the reconstructed sphere center coordinates and reconstructed sphere radius of the sphere in the single-layer sphere array calibration plate. The sphericity error is calculated based on the reconstructed sphere center coordinates and reconstructed sphere radius, and the radius error is calculated based on the reconstructed sphere radius. Based on the reconstructed sphere center coordinates, establish the sphere array coordinate direction matrix and calculate the sphere center array coordinates; Calculate the standard sphere center coordinates based on the standard sphere center coordinates in the physical coordinate system of the single-layer sphere array calibration plate; The array coordinate error is calculated based on the array coordinates of the sphere center and the standard sphere center-centered coordinates, and the sphere spacing error is calculated based on the standard sphere center-centered coordinates and the reconstructed sphere center coordinates. Define the ball count penalty and parameter regularization terms; The objective function value is obtained by weighted summing of the sphericity error, radius error, array coordinate error, sphere spacing error, sphere quantity penalty term, and parameter regularization term.

[0100] Furthermore, based on the three-dimensional reconstructed data, threshold segmentation and three-dimensional connected component analysis are performed to calculate the coordinates of the reconstructed sphere center and the radius of the reconstructed sphere in the single-layer sphere array calibration plate, including: The three-dimensional reconstructed volume data is segmented according to the set segmentation threshold to obtain binary volume data representing the spherical region; Based on the binary data, a three-dimensional connected component analysis was performed to obtain multiple spherical regions; Based on the voxel center coordinates of each voxel in the sphere region and the corresponding 3D reconstructed volume data, the coordinates of the reconstructed sphere center are calculated by gray-scale weighting. The radius of the reconstructed sphere is obtained by performing least-squares fitting on the boundary points of the sphere region.

[0101] Furthermore, the optimization module establishes a sphere array coordinate direction matrix based on the reconstructed sphere center coordinates, including: Obtain the coordinates of all reconstructed sphere centers to form a set of sphere centers; Calculate the sphere center covariance matrix based on the sphere center set, perform eigenvalue decomposition on the sphere center covariance matrix, and take the unit eigenvectors corresponding to the two largest eigenvalues ​​as the row and column directions of the sphere array coordinate system; The height direction is obtained by cross-multiplying the unit eigenvectors of the row and column directions. The row, column, and height directions together form the coordinate direction matrix of the spherical array.

[0102] Furthermore, the sphericity error is the sum of the squares of the normalized spherical residuals of the sphere, and the normalized spherical residual of a single sphere is the ratio of the root mean square residual of the distance from the boundary point of the sphere region to the center of the reconstructed sphere relative to the radius of the reconstructed sphere to the radius of the reconstructed sphere. The radius error is the mean square error between the reconstructed sphere radius and the standard sphere radius of all spheres; The sphere center array coordinates are the product of the transpose of the sphere array coordinate direction matrix and the reconstructed sphere center coordinates and the difference between the reconstructed sphere center coordinates and the mean of the reconstructed sphere center coordinates; The standard sphere center-centered coordinates are the difference between the standard sphere center coordinates and the mean of the standard sphere center coordinates in the physical coordinate system of the single-layer sphere array calibration plate. The array coordinate error is the mean square error of the difference between the array coordinates of the center of all spheres and the standard sphere center-centered coordinates; the sphere spacing error is the mean square error of the difference between the reconstructed sphere center spacing of all two adjacent spheres and the standard sphere center spacing. The ball count penalty term is the absolute value of the difference between the number of balls detected and the actual number of balls under the current candidate geometric parameters; the parameter regularization term is the sum of squared deviations of the changes between the initial geometric parameters and the current candidate parameters.

[0103] The CT system geometric parameter calibration method and apparatus based on reconstruction feedback optimization provided in the above embodiments have at least the following beneficial effects: (1) The optimization of geometric parameters is directly based on the three-dimensional reconstruction results, so that the calibration target no longer depends solely on the two-dimensional projection error, but rather the geometric parameters are evaluated jointly by the center position of the sphere array, the distance between the spheres, the sphericity and the radius consistency, thereby improving the accuracy, robustness and engineering applicability of CT geometric calibration; (2) The geometric parameters to be calibrated are clearly defined as the position of the X-ray source, the center position of the detector, the X-axis vector of the detector, and the Y-axis vector of the detector, so as to avoid treating the rotation angle, the geometry of each angle, the normal vector, the voxel projection point and the PCA principal direction as independent optimization parameters, thereby reducing the parameter dimension and parameter coupling. (3) Based on the CT system geometric coordinate system established on the rotation axis, it can clearly describe the spatial orientation of the X-ray source, detector center and detector pixel direction in the same coordinate system, so that the geometric model directly corresponds to the X-ray calculation process in CT reconstruction. (4) A single-layer M×N sphere array calibration plate is used. Geometric parameters are evaluated by the sphere center, sphere radius, sphericity and array spatial distribution in the reconstructed volume, so that the calibration target is directly oriented towards the final three-dimensional reconstruction quality, rather than relying solely on the two-dimensional projection domain error. (5) First, use PCA to construct the sphere array coordinate direction matrix from the reconstructed sphere center, then transform the reconstructed sphere center to the array coordinate system, and subtract it from the coordinates of the centered standard sphere array. This can reduce the influence of the calibration plate placement posture on the error evaluation. (6) The array coordinate error, sphere spacing error, sphericity error, radius error, sphere number penalty and parameter regularization term are combined into a unified objective function, so that the geometric parameter optimization takes into account the overall positional relationship, local sphere shape and equipment physical constraints, thus improving the stability of the calibration results.

[0104] Although preferred embodiments of the invention have been described, those skilled in the art, upon learning the basic inventive concept, can make other changes and modifications to these embodiments. Therefore, the appended claims are intended to be interpreted as including both the preferred embodiments and all changes and modifications falling within the scope of the invention. Clearly, those skilled in the art can make various alterations and modifications to the invention without departing from its spirit and scope. Thus, if these modifications and modifications of the invention fall within the scope of the claims and their equivalents, the invention is also intended to include these modifications and modifications.

Claims

1. A method for calibrating geometric parameters of a CT system based on reconstruction feedback optimization, characterized in that, include: Establish the geometric coordinate system of the CT system and determine the geometric parameters to be calibrated; Set up a single-layer ball array calibration board; The CT system is controlled to acquire projection images of a single-slice spherical array calibration plate at different projection angles according to the scanning trajectory, forming a set of projection images; The following operations are performed iteratively until a stopping condition is met: The position of each particle in the particle swarm is used as a set of candidate geometric parameters. Each projection image in the projection image set corresponds to a projection angle. Under the constraints of the geometric coordinate system of the CT system, 3D reconstruction is performed based on the projection image set and the corresponding candidate geometric parameters to obtain 3D reconstructed volume data corresponding to each set of candidate geometric parameters. Sphere data is extracted from the 3D reconstructed volume data. An objective function value with respect to error is solved based on the sphere data. The particle position is updated based on the individual historical optimal solution and the historical global optimal solution. Specifically, in the 3D reconstruction based on the projection image set and the corresponding candidate geometric parameters, the position of each particle is determined based on the candidate geometric parameters and the projection angle. The projection angle of each projected image in the image set is used to generate equivalent geometric parameters for each projected image: the three-dimensional position coordinates of the ray source, the three-dimensional position coordinates of the detector center, the three-dimensional displacement vector of the unit pixel in the X direction of the detector, and the three-dimensional displacement vector of the unit pixel in the Y direction of the detector are multiplied by the rotation matrix corresponding to the projection angle to obtain the equivalent geometric parameters of the projection angle; the equivalent geometric parameters of the projection angle include the three-dimensional position coordinates of the ray source, the three-dimensional position coordinates of the detector center, the three-dimensional displacement vector of the unit pixel in the X direction of the detector, and the three-dimensional displacement vector of the unit pixel in the Y direction of the detector under the projection angle. The globally optimal solution obtained when the stopping condition is met is used as the final calibration geometric parameters.

2. The method according to claim 1, characterized in that, In the geometric coordinate system of the CT system, the origin is the intersection of the plane containing the rotation axis and the center of rotation, the Z-axis is the rotation axis, and the X-axis and Y-axis are located in a plane perpendicular to the rotation axis and form a right-handed coordinate system with the Z-axis. In the geometric coordinate system of the CT system, the geometric parameters to be calibrated include the three-dimensional position coordinates of the X-ray source under the initial projection angle, the three-dimensional position coordinates of the detector center under the initial projection angle, the three-dimensional displacement vector per unit pixel in the X direction of the detector under the initial projection angle, and the three-dimensional displacement vector per unit pixel in the Y direction of the detector under the initial projection angle. A single-layer sphere array calibration board includes a substrate and spheres arranged in an array on the substrate.

3. The method according to claim 1, characterized in that, 3D reconstruction is performed based on the set of projected images and the corresponding candidate geometric parameters, including: The three-dimensional coordinates of the pixels in the projected image corresponding to the detector space are calculated based on the equivalent geometric parameters of the projection angle. The normal vector of the detector plane and the normal distance from the ray source to the detector plane are calculated based on the equivalent geometric parameters of the projection angle. The corresponding projected image is weighted by cone-beam averaging based on the pixel's three-dimensional coordinates and the normal distance from the X-ray source to the detector plane to obtain a weighted projected image. The weighted projection image is filtered along the X direction of the detector to obtain a filtered weighted projection image; Based on the normal vector under the projection angle and the normal distance from the X-ray source to the detector plane, the filtered weighted projection image is accumulated by voxel projection and back projection to obtain the three-dimensional reconstructed volume data corresponding to the candidate geometric parameters.

4. The method according to claim 3, characterized in that, Calculating the three-dimensional coordinates of pixels in the projected image corresponding to the detector space based on the equivalent geometric parameters of the projection angle includes: Calculate the center pixel coordinates of the detector based on its width and height; Calculate the offset of each pixel in the projected image relative to the coordinates of the detector's center pixel; Based on the offset, the three-dimensional displacement vector of a unit pixel in the X direction of the detector under the projection angle, and the three-dimensional displacement vector of a unit pixel in the Y direction of the detector under the projection angle, the displacement vector of the pixel in the detection space is calculated. Based on the three-dimensional position coordinates of the detector center under the projection angle and the displacement vector corresponding to the pixel in the detection space, the three-dimensional coordinates of each pixel in the projected image in the detector space are obtained. Calculate the normal vector under the projection angle and the normal distance from the ray source to the detector plane based on the equivalent geometric parameters of the projection angle, including: The normal vector under the projection angle is calculated based on the three-dimensional displacement vector of a unit pixel in the X direction of the detector under the projection angle and the three-dimensional displacement vector of a unit pixel in the Y direction of the detector under the projection angle. Based on the normal vector, the three-dimensional position coordinates of the X-ray source under the projection angle, and the three-dimensional position coordinates of the detector center under the projection angle, the normal distance from the X-ray source to the detector plane is calculated.

5. The method according to claim 3, characterized in that, Cone-beam weighting is performed on the corresponding projected image based on the pixel's three-dimensional coordinates and the normal distance from the X-ray source to the detector plane, including: The cone-beam weighting coefficients of the corresponding pixels are calculated based on the normal distance from the X-ray source to the detector plane and the three-dimensional coordinates of the pixels. The cone-beam weighting coefficients are applied to the corresponding pixels in the projected image to obtain a weighted projected image. Based on the normal vector at the projection angle and the normal distance from the X-ray source to the detector plane, the filtered weighted projection image is accumulated through voxel projection and back projection, including: Based on the normal vector, the three-dimensional position coordinates of the X-ray source under the projection angle, and the pre-determined center coordinates of the voxel, calculate the distance of the voxel relative to the X-ray source in the detector normal direction. If the distance of the voxel relative to the X-ray source in the detector normal direction is greater than 0, then the intersection point of the detector plane corresponding to the voxel is calculated based on the normal distance from the X-ray source to the detector plane, the distance of the voxel relative to the X-ray source in the detector normal direction, the three-dimensional position coordinates of the X-ray source under the projection angle, and the center coordinates of the voxel. The intersection points of the detector planes corresponding to the voxels are converted into continuous pixel coordinates on the filtered weighted projection image; Bilinear interpolation is performed at the continuous pixel coordinates on the filtered weighted projection image to obtain the back-projection image; The back-projected images under all projection angles are weighted and accumulated to obtain the 3D reconstructed volume data corresponding to the candidate geometric parameters.

6. The method according to claim 1, characterized in that, Extract sphere data from the 3D reconstructed volume data, and solve the objective function value with respect to the error based on the sphere data, including: Based on the three-dimensional reconstructed data, threshold segmentation and three-dimensional connected component analysis are performed to calculate the reconstructed sphere center coordinates and reconstructed sphere radius of the sphere in the single-layer sphere array calibration plate. The sphericity error is calculated based on the reconstructed sphere center coordinates and reconstructed sphere radius, and the radius error is calculated based on the reconstructed sphere radius. Based on the reconstructed sphere center coordinates, establish the sphere array coordinate direction matrix and calculate the sphere center array coordinates; Calculate the standard sphere center coordinates based on the standard sphere center coordinates in the physical coordinate system of the single-layer sphere array calibration plate; The array coordinate error is calculated based on the array coordinates of the sphere center and the standard sphere center-centered coordinates, and the sphere spacing error is calculated based on the standard sphere center-centered coordinates and the reconstructed sphere center coordinates. Define the ball count penalty and parameter regularization terms; The objective function value is obtained by weighted summing of the sphericity error, radius error, array coordinate error, sphere spacing error, sphere quantity penalty term, and parameter regularization term.

7. The method according to claim 6, characterized in that, Based on the 3D reconstructed data, threshold segmentation and 3D connected component analysis are performed to calculate the coordinates of the reconstructed sphere center and the radius of the reconstructed sphere in the single-layer sphere array calibration plate, including: The three-dimensional reconstructed volume data is segmented according to the set segmentation threshold to obtain binary volume data representing the spherical region; Based on the binary data, a three-dimensional connected component analysis was performed to obtain multiple spherical regions; Based on the voxel center coordinates of each voxel in the sphere region and the corresponding 3D reconstructed volume data, the coordinates of the reconstructed sphere center are calculated by gray-scale weighting. The radius of the reconstructed sphere is obtained by performing least-squares fitting on the boundary points of the sphere region.

8. The method according to claim 6, characterized in that, Based on the reconstructed sphere center coordinates, a sphere array coordinate direction matrix is ​​established, including: Obtain the coordinates of all reconstructed sphere centers to form a set of sphere centers; Calculate the sphere center covariance matrix based on the sphere center set, perform eigenvalue decomposition on the sphere center covariance matrix, and take the unit eigenvectors corresponding to the two largest eigenvalues ​​as the row and column directions of the sphere array coordinate system; The height direction is obtained by cross-multiplying the unit eigenvectors of the row and column directions. The row, column, and height directions together form the coordinate direction matrix of the spherical array.

9. The method according to claim 6, characterized in that, The sphericity error is the mean of the sum of squares of the normalized spherical residuals of the sphere, and the normalized spherical residual of a single sphere is the ratio of the root mean square residual of the distance from the boundary point of the sphere region to the center of the reconstructed sphere relative to the radius of the reconstructed sphere to the radius of the reconstructed sphere. The radius error is the mean square error between the reconstructed sphere radius and the standard sphere radius of all spheres; The sphere center array coordinates are the product of the transpose of the sphere array coordinate direction matrix and the reconstructed sphere center coordinates and the difference between the reconstructed sphere center coordinates and the mean of the reconstructed sphere center coordinates; The standard sphere center-centered coordinates are the difference between the standard sphere center coordinates and the mean of the standard sphere center coordinates in the physical coordinate system of the single-layer sphere array calibration plate. The array coordinate error is the mean square error of the difference between the array coordinates of the center of all spheres and the standard sphere center-centered coordinates; the sphere spacing error is the mean square error of the difference between the reconstructed sphere center spacing of all two adjacent spheres and the standard sphere center spacing. The ball count penalty term is the absolute value of the difference between the number of balls detected and the actual number of balls under the current candidate geometric parameters; the parameter regularization term is the sum of squared deviations of the changes between the initial geometric parameters and the current candidate parameters.

10. A calibrator for geometric parameters of a CT system based on reconstruction feedback optimization, applicable to any one of the methods described in claims 1-9, characterized in that, include: The coordinate establishment module is used to establish the geometric coordinate system of the CT system and determine the geometric parameters to be calibrated. The image acquisition module is used to control the CT system to acquire projected images of the single-layer spherical array calibration plate at different projection angles according to the scanning trajectory, forming a set of projected images; The optimization module iteratively performs the following operations until a stopping condition is met: taking the position of each particle in the particle swarm as a set of candidate geometric parameters, and under the constraints of the geometric coordinate system of the CT system, performing 3D reconstruction based on the set of projected images and the corresponding candidate geometric parameters to obtain 3D reconstructed volume data corresponding to each set of candidate geometric parameters; extracting sphere data based on the 3D reconstructed volume data, solving the objective function value with respect to the error based on the sphere data, and updating the particle positions; The parameter determination module is used to take the globally optimal solution obtained when the stopping condition is met as the final calibration geometric parameters.

Images