A high-precision industrial CT imaging system geometry parameter correction method and system

Abstract

The application discloses a high-precision industrial CT imaging system geometric parameter correction method and system, which comprises the following steps: step 1, a phantom is arranged at the center of a rotating table of an actual cone beam CT, the distance between a jth marker point in the phantom and a selected marker point is D real in a space coordinate system; step 2, actual projection data (v ij ,w ij ) of the jth marker point at an ith sampling angle in a detector local coordinate system is determined; step 3, predicted projection data (v i (ζ j ),w i (ζ j )) of the jth marker point at the ith sampling angle in the detector local coordinate system in a circular trajectory cone beam mode is obtained; step 4, the mean square error between (v ij ,w ij ) and (v i (ζ j ),w i (ζ j )) is set as a target function, and an optimal geometric parameter vector of the jth marker point is obtained, wherein the geometric parameter vector comprises the distance D rec between the jth marker point and the selected marker point in the space coordinate system and the distance sod between a ray source and the center of the rotating table; and step 5, the corrected distance sod real between the jth marker point and the selected marker point in the space coordinate system is obtained by using D rec , D real and sod.

Description

Technical Field

[0001] This invention relates to the field of CT (X-ray Computed Tomography) imaging technology, and in particular to a method and system for correcting the geometric parameters of a high-precision industrial CT imaging system. Background Technology

[0002] The essence of CT image artifacts lies in the discrepancy between the actual CT system and the ideal model. Artifacts caused by the deviation between the actual system's geometric parameters and the mathematical model parameters of the reconstruction algorithm are called geometric artifacts. Geometric artifacts often manifest in reconstructed images as blurred edges and misaligned object structures. Accurate CT system geometric parameters are indispensable for reconstructing high-quality images.

[0003] Due to mechanical errors in actual processing and installation, cone-beam CT systems may have problems such as difficulty in accurately locating the focal point of the X-ray source, the line connecting the focal point of the X-ray source and the rotation center of the examination bed not being perpendicular to the detector plane, and the detector placement not being consistent with the ideal position. Spiral cone-beam CT also needs to consider the problem that the forward direction of the examination bed is not parallel to the direction of rotation [2]. These problems will result in inaccurate geometric parameters of the system, which in turn will lead to geometric artifacts in the CT images.

[0004] Currently, there is a technique that uses a double-circular phantom to analytically solve for the geometric parameters of a CT system using only the projection data of the phantom at one angle. However, this method requires the introduction of some assumptions based on the double-circular phantom and places high demands on the placement and machining accuracy of the phantom.

[0005] Another existing method uses a phantom composed of four metal spheres. Given the distances between adjacent spheres within the phantom, a least-squares objective function is constructed using the projected addresses and predicted coordinates of the spheres. Then, optimization methods are used to solve for the geometric parameters of the CT system. While this method can be applied to circular, spiral, or other scanning trajectories, it requires continuous image reconstruction, making the calculations quite complex. Summary of the Invention

[0006] The purpose of this invention is to provide a method and system for correcting the geometric parameters of a high-precision industrial CT imaging system, so as to overcome or at least mitigate at least one of the above-mentioned defects of the prior art.

[0007] To achieve the above objectives, the present invention provides a method for correcting the geometric parameters of a high-precision industrial CT imaging system, characterized by comprising:

[0008] Step 1: Place the phantom at the center of the turntable of the actual cone-beam CT. The phantom has one marker point, where the distance between the j-th marker point and a selected marker point in the spatial coordinate system is D. real The spatial coordinate system has its origin at the center of the turntable and its axis vectors are the same as those of the world coordinate system.

[0009] Step 2, determine the actual projection data (v) of the j-th marker point at the i-th sampling angle in the detector's local coordinate system. ij ,w ij The detector's local coordinate system has its origin at the detector center, and its three coordinate axes are rotated clockwise by γ and θ around the x, y, and z axes, respectively. y ,θ z Spend;

[0010] Step 3, obtain the predicted projection data (v) of the j-th marker point in the detector's local coordinate system at the i-th sampling angle under the circular trajectory cone-beam configuration. i (ζ j ),w i (ζ j ));

[0011] Step 4, transfer the actual projection data (v) from Step 2. ij ,w ij ) and the predicted projection data (v) from step 3 i (ζ j ),w i (ζ j Mean square error between )) The objective function is set as fun, and the optimal geometric parameter vector of the j-th marker point is obtained using an optimization method. The geometric parameter vector includes the distance D between the j-th marker point and the selected marker point in the spatial coordinate system. rec And the distance between the radiation source and the center of the turntable (sod);

[0012] Step 5, use D from step 1 real and step 4, D rec And sod, obtain the corrected distance sod between the j-th marker point and the selected marker point in the spatial coordinate system. real .

[0013] Furthermore, the optimization method takes minimizing fun as the optimization objective, as shown in equation (1):

[0014]

[0015] In the formula, ζ j Let ζ represent the geometric parameter vector of the j-th marker point. up ,ζ down Representing ζ respectively j The upper and lower bounds of F(ζ)up ,ζ down ) represents the penalty function term, and K represents the number of sampling angles.

[0016] Furthermore, the optimization method in step 4 specifically includes:

[0017] Step 41, convert the n*1 dimensional geometric parameter vector ζ j The initial value vector is set to X0; where n represents ζ. j The number of parameters included;

[0018] Step 42: Calculate the objective function fun(X0) at X0, assign X0 to the current geometric parameter vector current_sol and the historical best geometric parameter vector best_sol, and assign fun(x0) to fun(current_sol) and fun(best_sol), where fun(current_sol) = current_cost and fun(best_sol) = best_cost.

[0019] Step 43: Use the current_sol from step 42 as the initial value vector, and use the first simplex method to obtain the new solution current_sol and current_cost. If fun(current_sol) is less than best_cost, then assign the value of current_sol to best_sol and the value of fun(current_sol) to best_cost.

[0020] Step 44: At current_sol, use the improved neighborhood search method to search for a random point in a certain neighborhood as a new initial value vector new_sol: If fun(new_sol) is less than current_cost, then assign new_sol to current_sol; if fun(new_sol) is not less than current_cost, and the current loop state has reached the preset loop requirement, then assign new_sol to current_sol; if fun(new_sol) is not less than current_cost, and the current loop state has not reached the preset loop requirement, then execute step 43 until best_sol is less than the first preset threshold.

[0021] Furthermore, step 43, the first simplex method, specifically includes:

[0022] Step 431: Based on the initial value vector, create n linearly independent vertices X. p =X0+Q p *e p Q pe represents the step size corresponding to the p-th vertex. p Let current_sol represent the n-dimensional unit vector corresponding to the p-th vertex. If current_sol is less than the second preset threshold that can be obtained through multiple iterations of the simplex method, then the step size of the key parameters is restricted to make convergence faster.

[0023] Step 432: When the optimization method loops to the preset number of iterations of the first simplex method, the second simplex method is used intermittently to find a solution that meets the requirements as the initial value vector for the next iteration of the first simplex method; wherein, the preset number of iterations of the second simplex method is less than that of the first simplex method.

[0024] Furthermore, the marker point is the pre-set center of mass corresponding to the ball.

[0025] This invention also provides a high-precision calibration system for geometric parameters of an industrial CT imaging system, comprising:

[0026] The phantom, set at the center of the turntable of the actual cone-beam CT, has l marker points, where the distance between the j-th marker point and a selected marker point in the spatial coordinate system is D. real The spatial coordinate system has its origin at the center of the turntable and its axis vectors are the same as those of the world coordinate system.

[0027] The actual projection data acquisition unit is used to determine the actual projection data (v) of the j-th marker point at the i-th sampling angle in the detector's local coordinate system. ij ,w ij The detector's local coordinate system has its origin at the detector center, and its three coordinate axes are rotated clockwise by γ and θ around the x, y, and z axes, respectively. y ,θ z Spend;;

[0028] The predictive projection data acquisition unit is used to acquire the predictive projection data (v) of the j-th marker point in the detector's local coordinate system at the i-th sampling angle under the circular trajectory cone-beam configuration. i (ζ j ),w i (ζ j ));

[0029] The geometric parameter optimization unit is used to obtain the actual projection data (v) from the actual projection data acquisition unit. ij ,w ij ) and the predicted projection data (v) of the predicted projection data acquisition unit i (ζ j ),w i (ζ j Mean square error between )) The objective function is set as fun, and the optimal geometric parameter vector of the j-th marker point is obtained using an optimization method. The geometric parameter vector includes the distance D between the j-th marker point and the selected marker point in the spatial coordinate system. rec And the distance between the radiation source and the center of the turntable (sod);

[0030] Distance correction unit, which is used for D of the phantom real and the D of the geometric parameter optimization unit rec And sod, obtain the corrected distance sod between the j-th marker point and the selected marker point in the spatial coordinate system. real .

[0031] Furthermore, the optimization method takes minimizing fun as the optimization objective, as shown in equation (1):

[0032]

[0033] In the formula, ζ j Let ζ represent the geometric parameter vector of the j-th marker point. up ,ζ down Representing ζ respectively j The upper and lower bounds of F(ζ) up ,ζ down ) represents the penalty function term, and K represents the number of sampling angles.

[0034] Furthermore, the optimization method for the geometric parameter optimization unit specifically includes:

[0035] The initial value vector is set in the sub-unit, which is used to set the n*1 dimensional geometric parameter vector ζ. j The initial value vector is set to X0; where n represents ζ. j The number of parameters included;

[0036] The initial value calculation sub-unit is used to calculate the objective function fun(X0) at X0. X0 is assigned to the current geometric parameter vector current_sol and the historical best geometric parameter vector best_sol. fun(x0) is assigned to fun(current_sol) and fun(best_sol). fun(current_sol) = current_cost, fun(best_sol) = best_cost.

[0037] The simplex method sub-element is used to take the current_sol of the initial value calculation sub-element as the initial value vector, and use the first simplex method to obtain the new solution current_sol and current_cost. If fun(current_sol) is less than best_cost, then current_sol is assigned to best_sol and fun(current_sol) is assigned to best_cost.

[0038] The neighborhood search subunit is used to search for a random point in a neighborhood at current_sol using an improved neighborhood search method as a new initial value vector new_sol: if fun(new_sol) is less than current_cost, then new_sol is assigned to current_sol; if fun(new_sol) is not less than current_cost, and the current loop state has reached the preset loop requirement, then new_sol is assigned to current_sol; if fun(new_sol) is not less than current_cost, and the current loop state has not reached the preset loop requirement, then the process continues until the best_sol obtained by the simplex method subunit is less than the first preset threshold.

[0039] Furthermore, the simplex subunit specifically includes:

[0040] The vertex creation module is used to create n linearly independent vertices X based on the initial value vector. p =X0+Q p *e p Q p e represents the step size corresponding to the p-th vertex. p Let current_sol represent the n-dimensional unit vector corresponding to the p-th vertex. If current_sol is less than the second preset threshold that can be obtained through multiple iterations of the simplex method, then the step size of the key parameters is restricted to make convergence faster.

[0041] The loop count determination module is used to intermittently use the second simplex method when the optimization method reaches the preset loop count of the first simplex method in order to find a solution that meets the requirements as the initial value vector for the next iteration of the first simplex method; wherein the preset loop count of the second simplex method is less than that of the first simplex method.

[0042] Furthermore, the marker point is the pre-set center of mass corresponding to the ball.

[0043] The present invention has the following advantages due to the adoption of the above technical solutions:

[0044] 1. This invention uses the Nelder-Mead method repeatedly, that is, after the Nelder-Mead method is used continuously, a new solution is randomly searched in its neighborhood as the initial point for the next Nelder-Mead method. Therefore, this invention can overcome the problem of getting trapped in local minima due to poor initial value vector.

[0045] 2. This invention uses the alternating search of large and small simplexes and employs a penalty function to keep the simplex within a reasonable range. Therefore, this invention can solve the problems of excessively long search times and instability when looking for solutions that meet expectations.

[0046] 3. After finding a relatively good solution, this invention restricts the search range of certain parameters that have a significant impact on the objective function, which can lead to faster convergence. Furthermore, when creating the vertices of the simplex, the change of these parameters is limited. Therefore, this invention can obtain a better solution and converge faster. Attached Figure Description

[0047] Figure 1 This is a schematic diagram of the mold according to an embodiment of the present invention. Detailed Implementation

[0048] In the accompanying drawings, the same or similar reference numerals are used to denote the same or similar elements or elements having the same or similar functions. The embodiments of the present invention will now be described in detail with reference to the accompanying drawings.

[0049] In the description of this invention, the terms "center," "longitudinal," "lateral," "front," "rear," "left," "right," "vertical," "horizontal," "top," "bottom," "inner," and "outer," etc., indicate the orientation or positional relationship based on the orientation or positional relationship shown in the accompanying drawings. They are used only for the convenience of describing this invention and simplifying the description, and do not indicate or imply that the device or element referred to must have a specific orientation, or be constructed and operated in a specific orientation. Therefore, they should not be construed as limiting the scope of protection of this invention.

[0050] The method for correcting the geometric parameters of a high-precision industrial CT imaging system provided in this invention includes:

[0051] Step 1: Set the phantom at the center of the turntable of the actual cone-beam CT. The phantom has l marker points, where the distance between the j-th marker point and a selected marker point in the spatial coordinate system is D. real .like Figure 1 As shown, this embodiment uses a model composed of two forest spheres (rubies) at a known distance, i.e., l = 2. The marker point is the pre-set centroid corresponding to the small sphere, that is, D real Given parameters, and the distance between the centroids of the two forest spheres.

[0052] In this system, the spatial coordinate system has its origin at the center of the turntable and shares the same axis vectors as the world coordinate system. The world coordinate system axes include the x, y, and z axes. The x-axis is the axis of the line connecting the radiation source and the detector; the z-axis is the axis perpendicular to the horizontal plane and pointing upwards; the y-axis is perpendicular to both the x and z axes. Rotation around the x-axis yields γ, and rotation around the y-axis yields θ. y θ is obtained by rotating around the z-axis. z .

[0053] Step 2, determine the actual projection data (v) of the j-th marker point at the i-th sampling angle in the detector's local coordinate system. ij ,w ij The detector's local coordinate system has its origin at the detector center, and its three coordinate axes are rotated clockwise by γ and θ around the x, y, and z axes, respectively. y ,θ z Spend.

[0054] Step 3, obtain the predicted projection data (v) of the j-th marker point in the detector's local coordinate system at the i-th sampling angle under the circular trajectory cone-beam configuration. i (ζ j ),w i (ζ j )).

[0055] Step 4, transfer the actual projection data (v) from Step 2. ij ,w ij ) and the predicted projection data (v) from step 3 i (ζ j ),w i (ζ j Mean square error between )) The objective function is set as fun, and the optimal geometric parameter vector of the j-th marker point is obtained using an optimization method. The geometric parameter vector includes the distance D between the j-th marker point and the selected marker point in the spatial coordinate system. rec The distance between the X-ray source and the center of the turntable (sod); the distance between the intersection point of the straight line passing through the center of the turntable and the detector and the X-ray source (sdd); and the offset of the detector center from its ideal position within the detector plane (uv). u ,offset v ), the rotation angle of the detector relative to the world coordinate system (θ) y ,θ z ,γ), and the coordinates (x, y) of the j-th marker point in the spatial coordinate system. j ,y j ,z j ).

[0056] Step 5, use D from step 1 real and step 4, D rec And SOD, calculated using the formula: Dreal / sod rec =D rec / sod, obtain the corrected distance sod between the j-th marker point and the selected marker point in the spatial coordinate system. rec The optimal geometrical vector of the corrected j-th marker point can be used for reconstruction to achieve the desired application.

[0057] In one embodiment, the optimization method aims to minimize fun, as shown in equation (1):

[0058]

[0059] In the formula, ζ j Let ζ represent the geometric parameter vector of the j-th marker point. up ,ζ down Representing ζ respectively j The upper and lower bounds of F(ζ) up ,ζ down ,ζ j ) represents the penalty function term. For example, it is equal to 0 when it does not exceed the upper or lower bound, and is assigned a value of approximately infinity when it exceeds the upper or lower bound. K represents the number of sampling angles. In this method, setting K=60 has already achieved good results.

[0060] Of course, optimization methods can also achieve other optimization objectives besides minimization, depending on actual needs.

[0061] In one embodiment, the optimization method in step 4 specifically includes:

[0062] Step 41, convert the n*1 dimensional geometric parameter vector ζ j The initial value vector is set to X0. Here, n represents ζ. j The number of parameters included, such as ζ j When the vector is a 1*13 row vector, n is 13.

[0063] Step 42: Calculate the objective function fun(X0) at X0. Assign X0 to the current geometric parameter vector current_sol and the historical best geometric parameter vector best_sol. Assign fun(x0) to fun(current_sol) and fun(best_sol). fun(current_sol) = current_cost, fun(best_sol) = best_cost.

[0064] Step 43: Using current_sol from step 42 as the initial value vector, obtain the new solution current_sol and current_cost using the first simplex method. If fun(current_sol) is less than best_cost, then assign the value of current_sol to best_sol and the value of fun(current_sol) to best_cost.

[0065] Step 44: At current_sol, use the improved neighborhood search method to search for a random point in a certain neighborhood as a new initial value vector new_sol: If fun(new_sol) is less than current_cost, then assign new_sol to current_sol; if fun(new_sol) is not less than current_cost, and the current loop state has reached the preset loop requirement, such as executing the operation every n times, then assign new_sol to current_sol; if fun(new_sol) is not less than current_cost, and the current loop state has not reached the preset loop requirement, then execute step 43 until best_sol is less than the first preset threshold, such as about 1 / 5 of the detector unit length.

[0066] The improved neighborhood search method in this embodiment has a strong ability to escape local minima before a better solution is found. When a better solution is found (i.e., when best_sol is relatively small), it can limit the random search range of certain key geometric parameters to make the objective function value of the updated solution smaller, so that a solution with a smaller objective function can be obtained after using the simplex method.

[0067] In one embodiment, the method of "using the current_sol from step 42 as the initial value vector to obtain the new solution current_sol" in step 43 specifically includes:

[0068] Step 431: Based on the initial value vector, create n linearly independent vertices X. p =X0+Q p *e p Q p It is a number randomly selected from a certain range, representing the step size corresponding to the p-th vertex, e. pLet `current_sol` represent the n-dimensional unit vector corresponding to the p-th vertex. The subsequent steps are consistent with the standard simplex method. If `current_sol` is less than a second preset threshold that can be obtained through multiple simplex iterations (e.g., a pixel size), then limiting the step size of the key parameters (such as `offset_x` and `γ`) makes convergence faster. This can, to some extent, replace random neighborhood search to escape local minima. Furthermore, this method allows for the use of better points obtained in the previous simplex iteration as initial value vectors, resulting in better efficiency.

[0069] Step 432: When the optimization method iterates to the preset number of iterations of the first simplex method, the second simplex method is interspersed. For example, the second simplex method can be used every n iterations. The purpose of setting the second simplex method to iterate here is to replace the random neighborhood search to a certain extent, so as to escape the search for local minima, and to find a solution with a smaller objective function value as the initial value vector for the next iteration of the first simplex method. The preset number of iterations for the second simplex method is less than that for the first simplex method. In this embodiment, a larger random step size can be set at the second simplex method to escape local minima, and a shorter step size can be set before some parameters that have a significant impact on the objective function when the initial value vector of the first simplex method meets the requirements, as shown in step 431 above. The second simplex method can balance speed and obtain a better solution than random neighborhood search.

[0070] In this embodiment, simplexes of random sizes are used continuously. This balances the ability to quickly find a better solution and escape local optima to a certain extent. It also incorporates a flexible neighborhood random search algorithm. When current_sol is less than one pixel, the solution is considered good enough. At this point, some parameters that have a significant impact on the objective function value can be restricted.

[0071] This invention also provides a high-precision geometric parameter correction system for an industrial CT imaging system, comprising a phantom, an actual projection data acquisition unit, a predicted projection data acquisition unit, a geometric parameter optimization unit, and a distance correction unit, wherein:

[0072] The phantom is positioned at the center of the turntable of an actual cone-beam CT scanner. The phantom has one marker point, where the distance between the j-th marker point and a selected marker point in the spatial coordinate system is D. real The spatial coordinate system has its origin at the center of the turntable and its axis vectors are the same as those of the world coordinate system.

[0073] The actual projection data acquisition unit is used to determine the actual projection data (v) of the j-th marker point at the i-th sampling angle in the detector's local coordinate system. ij ,w ijThe detector's local coordinate system has its origin at the detector center, and its three coordinate axes are rotated clockwise by γ and θ around the x, y, and z axes, respectively. y ,θ z Spend.

[0074] The predictive projection data acquisition unit is used to acquire the predictive projection data (v) of the j-th marker point in the detector's local coordinate system at the i-th sampling angle under the circular trajectory cone-beam configuration. i (ζ j ),w i (ζ j )).

[0075] The geometric parameter optimization unit is used to obtain the actual projection data (v) from the actual projection data acquisition unit. ij ,w ij ) and the predicted projection data (v) of the predicted projection data acquisition unit i (ζ j ),w i (ζ j Mean square error between )) The objective function is set as fun, and the optimal geometric parameter vector of the j-th marker point is obtained using an optimization method. The geometric parameter vector includes the distance D between the j-th marker point and the selected marker point in the spatial coordinate system. rec And the distance between the radiation source and the center of the turntable, sod.

[0076] The distance correction unit is used for the D of the phantom real and the D of the geometric parameter optimization unit rec And sod, obtain the corrected distance sod between the j-th marker point and the selected marker point in the spatial coordinate system. real .

[0077] In one embodiment, the optimization method of the geometric parameter optimization unit specifically includes an initial value vector setting subunit, an initial value calculation subunit, a simplex method subunit, and a neighborhood search subunit, wherein:

[0078] The initial value vector setting sub-unit is used to set the n*1 dimensional geometric parameter vector ζ. j The initial value vector is set to X0; where n represents ζ. j The number of parameters included.

[0079] The initial value calculation subunit is used to calculate the objective function fun(X0) at X0. X0 is assigned to the current geometric parameter vector current_sol and the historical best geometric parameter vector best_sol. fun(x0) is assigned to fun(current_sol) and fun(best_sol). fun(current_sol) = current_cost, fun(best_sol) = best_cost.

[0080] The simplex method sub-element is used to take the current_sol of the initial value calculation sub-element as the initial value vector, and use the first simplex method to obtain the new solution current_sol and current_cost. If fun(current_sol) is less than best_cost, then current_sol is assigned to best_sol and fun(current_sol) is assigned to best_cost.

[0081] The neighborhood search subunit is used to search for a random point in a neighborhood at current_sol using the improved neighborhood search method as a new initial value vector new_sol: if fun(new_sol) is less than current_cost, then new_sol is assigned to current_sol; if fun(new_sol) is not less than current_cost, and the current loop state has reached the preset loop requirement, then new_sol is assigned to current_sol; if fun(new_sol) is not less than current_cost, and the current loop state has not reached the preset loop requirement, then the process continues until the best_sol obtained by the simplex method subunit is less than the first preset threshold.

[0082] In one embodiment, the simplex sub-unit specifically includes a vertex creation module and a loop count determination module, wherein:

[0083] The vertex creation module is used to create n linearly independent vertices X based on the initial value vector. p =X0+Q p *e p Q p This represents the step size corresponding to the p-th vertex, where p can be set to a random number within a reasonable range, and e p Let current_sol represent the n-dimensional unit vector corresponding to the p-th vertex. If current_sol is less than the second preset threshold that can be obtained through multiple iterations of the simplex method, then limiting the step size of the key parameters will make convergence faster.

[0084] The iteration count determination module is used to intermittently use the second simplex method when the optimization method iterates to the preset iteration count of the first simplex method in order to find a solution that meets the requirements as the initial value vector for the next iteration of the first simplex method; wherein, the preset iteration count of the second simplex method is less than that of the first simplex method.

[0085] The advantages of this invention are: First, it can find a satisfactory solution in a very short time and with sufficient stability. Second, it can minimize the dependence on the initial value vector after determining an acceptable range of penalty functions. Furthermore, this method not only requires readily available and simple modalities, but also allows for the recovery of a relatively accurate solution with minimal sampling angles.

[0086] The feasibility of this method is verified through simulation experiments. In the simulated circular trajectory scanning, the detector resolution is 1536*1536, the detector element size is 0.278mm, the number of samples is 60, and the sampling interval is 6°. The phantom used in this simulation experiment is as shown above: two ruby ​​spheres with a radius of 2mm, with a horizontal distance of 25mm and a vertical distance of 20mm. During orthographic projection, the number of voxels is 1536*1536*1351, and the image size is 46.1612*46.1612*40.6829. The geometric parameters of distance in the table below are in mm, and the geometric parameters of angle are in degrees. The penalty function range is:

[0087]

[0088]

[0089] The initial value vector is:

[0090] The objective function represents the root mean square error between the actual and predicted projections of the 60 centroids, which can be considered as the Euclidean distance between them. In this case, the objective function is 0.0560. The distance error is approximately 0.2 pixels.

[0091] The size of the prime element can be considered an acceptable solution.

[0092]

[0093] Under the condition of being consistent with the above true values, one thousand initial value vectors were randomly selected from the intervals sdd:[1555,1645], sod:[155,245], x00 / y00 / x10 / y10 / z00 / z10:[-45,45], offset_u / offset_v:[-3.6,3.6], theta_y\theta_z\gamma:[-1.8,1.8] which follow a normal distribution. The simplex and large loop were terminated when the threshold was set to 0.0561. The total time was 1906 seconds, and the average time to find a solution was 1.9 seconds each time.

[0094] This invention proposes an optimized algorithm for further searching the global minimum solution based on the objective function used in the present method, namely, the improved simplex method. This improved simplex method improves the ability and speed of searching for the global optimum by continuously using simplex methods with different iteration numbers and sizes, and by setting a more flexible penalty function. Furthermore, this invention incorporates a small number of restricted neighborhood random searches during the iteration process, depending on the situation. This allows for more effective escape from the minimum solution while receiving the best possible solution, so that it can be used as the initial value vector for the next simplex method, resulting in better performance.

[0095] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and not to limit them. Those skilled in the art should understand that modifications can be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features; these modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method for calibrating the geometric parameters of a high-precision industrial CT imaging system, characterized in that, Comprising: Step 1, set a phantom in the center of the gantry of the actual cone beam CT, the phantom has 1 marker point, wherein the distance between the jth marker point and a selected marker point in the spatial coordinate system is D real ; wherein, the spatial coordinate system takes the center of the gantry as the origin and has the same axial vector as the world coordinate system; Step 2, determine the actual projection data of the jth marker point at the ith sampling angle in the detector local coordinate system (v ij , ij ); wherein the detector local coordinate system takes the detector center as the origin, and its three coordinate axes are rotated clockwise by γ, θ y , z degrees around the x, y, z axes respectively; Step 3, obtaining the predicted projection data of the jth marker point at the ith sampling angle in the detector local coordinate system under the circular trajectory C-arm mode (v i (ζ j ),w i (ζ j )) in the detector local coordinate system, where ζ j represents the geometric parameter vector of the jth marker point. Step 4, transfer the actual projection data (v) from Step 2. ij ,w ij ) and the predicted projection data (v) from step 3 i (ζ j ),w i (ζ j Mean square error between )) The objective function is set as fun, and the optimal geometric parameter vector of the j-th marker point is obtained using an optimization method. The geometric parameter vector includes the distance D between the j-th marker point and the selected marker point in the spatial coordinate system. rec And the distance between the radiation source and the center of the turntable (sod); Step 5, using D of Step 1 real , and D of Step 4 rec and sod, to obtain the corrected distance sod of the jth marker point and the selected marker point in the spatial coordinate system real ; The optimization method of step 4 specifically comprises: Step 41, set the initial vector of the n*l dimensional geometric parameter vector ζ j to X0; wherein n represents the number of parameters contained in ζ j ; Step 42, calculating the objective function fun(X0) at X0, assigning X0 to the current geometric parameter vector current_sol and the historically optimal geometric parameter vector best_sol, assigning fun(x0) to fun(current_sol) and fun(best_sol), fun(current_sol) = current_cost, and fun(best_sol) = best_cost; Step 43, taking the current_sol of step 42 as the initial value vector, using the first simplex method to obtain a new solution current_sol and current_cost, and if fun(current_sol) is less than best_cost, assigning current_sol to best_sol and assigning fun(current_sol) to best_cost; Step 44, using the improved neighborhood search method at current_sol to search for a random point in a certain neighborhood as a new initial value vector new_sol: if fun(new_sol) is less than current_cost, assigning new_sol to current_sol; if fun(new_sol) is not less than current_cost and the current loop state has reached the preset loop requirement, assigning new_sol to current_sol; if fun(new_sol) is not less than current_cost and the current loop state has not reached the preset loop requirement, executing step 43 until best_sol is less than the first preset threshold.

2. The method of claim 1, wherein the method is performed by a computer system. The optimization method takes minimizing fun as the optimization objective, as shown in equation (1): In the formula, ζ up ,ζ down respectively represent the upper and lower bounds of ζ j , F(ζ up ,ζ down ) represents a penalty function term, and K represents the number of sampling angles.

3. The method of claim 1, wherein the method is performed by a computer system. The first simplex method of step 43 specifically comprises: Step 431, according to the initial value vector, create n linearly independent vertex X p = X0+ Q p *e p , Q p indicates the step size corresponding to the pth vertex, e p indicates the n-dimensional unit vector corresponding to the pth vertex, if current_sol is less than the second preset threshold value that can be obtained through multiple simplex method iterations, the step size of the key parameter is limited to make the convergence faster; Step 432, when the optimization method loops to the preset loop number of the first simplex method, the second simplex method is used to find a solution that meets the requirements as the initial value vector of the next first simplex method; wherein the preset loop number of the second simplex method is less than that of the first simplex method.

4. The method of claim 1, wherein the method is performed by a computer system. The marker point is the centroid corresponding to the pre-set small ball.

5. A high-precision correction system for geometric parameters of an industrial CT imaging system, characterized in that, Comprising: A phantom is arranged at the center of a rotation table of an actual cone beam CT, and the phantom has 1 marker point, wherein the distance between the jth marker point and a selected marker point in a spatial coordinate system is D real ; wherein the spatial coordinate system takes the center of the rotation table as an origin and has the same axial vector as the world coordinate system. an actual projection data acquisition unit, configured to determine actual projection data (v ij , ij ) of the jth marker point at the ith sampling angle in a detector local coordinate system; wherein the detector local coordinate system takes the center of the detector as the origin, and three coordinate axes thereof are rotated clockwise by γ, θ y , z degrees about the x, y, and z axes, respectively. a prediction projection data acquisition unit, configured to acquire predicted projection data of a jth marker point at an ith sampling angle in a local coordinate system of a detector in a circular trajectory cone beam mode (v i (ζ j ),w i (ζ j )) and ζ j denotes a geometric parameter vector of the jth marker point; a geometry parameter optimization unit configured to set a mean square error between actual projection data (v ij , w ij ) of the actual projection data acquisition unit and predicted projection data (v i (ζ j , w i (ζ j )) of the predicted projection data acquisition unit as a target function fun, and obtain an optimal geometry parameter vector of the jth marker point by using an optimization method, the geometry parameter vector including a distance D rec of the jth marker point to the selected marker point in a spatial coordinate system and a distance sod of the ray source to the center of the turntable; a geometry parameter optimization unit configured to set a mean square error between actual projection data (v ij , w ij ) of the actual projection data acquisition unit and predicted projection data (v i (ζ j , w i (ζ j )) of the predicted projection data acquisition unit as a target function fun, and obtain an optimal geometry parameter vector of the jth marker point by using an optimization method, the geometry parameter vector including a distance D rec of the jth marker point to the selected marker point in a spatial coordinate system and a distance sod of the ray source to the center of the turntable; a distance correction unit configured to correct the distance between the jth marker point and the selected marker point in the spatial coordinate system based on the D real of the phantom and the D rec of the geometry parameter optimization unit, to obtain a corrected distance sod real between the jth marker point and the selected marker point in the spatial coordinate system The optimization method of the geometric parameter optimization unit specifically comprises: a initial value vector setting sub-unit, configured to set a initial value vector of the n*1 dimensional geometric parameter vector ζ j as X0; wherein, n represents the number of parameters contained in ζ j ; The initial value calculation sub-unit is configured to calculate the objective function fun(X0) at X0, assign X0 to the current geometric parameter vector current_sol and the historically optimal geometric parameter vector best_sol, assign fun(x0) to fun(current_sol) and fun(best_sol), fun(current_sol) = current_cost, and fun(best_sol) = best_cost; The simplex method subunit is configured to use the current_sol of the initial value calculation subunit as an initial value vector to obtain a new solution current_sol and a current_cost by using a first simplex method, and if the fun(current_sol) is less than the best_cost, assign the current_sol to the best_sol and assign the fun(current_sol) to the best_cost. The neighborhood search subunit is configured to search for a random point as a new initial value vector new_sol in a neighborhood at the current_sol by using an improved neighborhood search method, and if the fun(new_sol) is less than the current_cost, assign the new_sol to the current_sol, if the fun(new_sol) is not less than the current_cost and the current loop state has reached a preset loop requirement, assign the new_sol to the current_sol, and if the fun(new_sol) is not less than the current_cost and the current loop state has not reached the preset loop requirement, until the best_sol obtained by the simplex method subunit is less than a first preset threshold.

6. The high-precision industrial CT imaging system geometry parameter correction system of claim 5, wherein, The optimization method takes minimizing the fun as an optimization objective, as shown in equation (1): In the formula, ζ up ,ζ down respectively represent the upper and lower bounds of ζ j , F(ζ up ,ζ down ) represents a penalty function term, and K represents the number of sampling angles.

7. The high-precision industrial CT imaging system geometry parameter correction system of claim 5, wherein, The simplex method subunit specifically includes: a vertex creation module, configured to create n linearly independent vertices X according to the initial value vector p = X0+ Q p e p , Q p denotes a step size corresponding to the pth vertex, e p denotes an n-dimensional unit vector corresponding to the pth vertex, and if the current_sol is less than a second preset threshold value that can be obtained through multiple simplex method iterations, the step size of the key parameter is limited to make the convergence faster The loop number judgment module is configured to, when the optimization method loops to a preset loop number of the first simplex method, use a second simplex method to find a solution meeting a requirement as an initial value vector of the next first simplex method, and the preset loop number of the second simplex method is less than that of the first simplex method.

8. The high-precision industrial CT imaging system geometry parameter correction system of claim 5, wherein, The marker point is a centroid corresponding to a pre-set small ball.

Images