Method and system for correcting geometric parameters of high-precision industrial CT imaging system

Abstract

The invention discloses a high-precision industrial CT imaging system geometric parameter correction method and system, and the method comprises the steps: 1, enabling a die body to be disposed at the center of a rotary table of an actual cone beam CT, and enabling the distance between a jth mark point in the die body and a selected mark point in a space coordinate system to be Dreal; step 2, determining actual projection data (vij, wij) of the jth mark point at the ith sampling angle under the local coordinate system of the detector; step 3, acquiring predicted projection data (vi (zeta j), wi (zeta j)) of a jth mark point at an ith sampling angle under a local coordinate system of the detector in a circular track cone beam mode; 4, a mean square error between (vij, wij) and (vi (zeta j), wi (zeta j)) is set as a target function, the optimal geometric parameter vector of the jth mark point is obtained, and the geometric parameter vector comprises the distance Drec between the jth mark point and the selected mark point in the space coordinate system and the distance sod between the ray source and the center of the rotary table; and step 5, using the Dreal, the Drec and the sod to obtain the correction distance sreal of the jth mark point and the selected mark point in the space coordinate system.

Description

Technical Field

[0001] The present invention relates to the field of CT (full name in English is "X-ray Computed Tsomography", full name in Chinese is "X-ray computer tomography") imaging technology, and in particular to a high-precision correction method and system for geometric parameters of an industrial CT imaging system. Background Art

[0002] The essence of CT image artifacts is the difference between the actual CT system and the ideal model. The artifacts caused by the deviation between the actual system geometric parameters and the reconstruction algorithm mathematical model parameters are called geometric artifacts. Geometric artifacts often appear in the reconstructed image 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, the cone-beam CT system may have problems such as the difficulty in accurately locating the focal point of the radiation source, the line connecting the focal point of the radiation source and the rotation center of the examination bed is not perpendicular to the detector plane, and the detector placement is inconsistent 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 rotation axis [2]. These problems will lead to inaccurate system geometric parameters, which will in turn cause geometric artifacts in the CT image.

[0004] Currently, there is an existing technology that can use a double-circular phantom to analyze and solve the geometric parameters of the CT system with only the projection data of the phantom at one angle. However, this method requires the introduction of some assumptions with the help of the double-circular phantom, and has high requirements for the placement and processing accuracy of the phantom.

[0005] Another existing method is to use a phantom composed of four metal balls. Under the premise of knowing the distance between adjacent balls in the phantom, the least squares objective function is constructed by the collected projection addresses of the balls and the predicted coordinates of the balls, and then the geometric parameters of the CT system are solved using the optimization method. Although this method can be applied to circular trajectories, spiral trajectories or other scanning trajectories, it requires continuous image reconstruction and the calculation is relatively complex. Summary of the invention

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

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

[0008] Step 1: Set the phantom at the center of the turntable of the actual cone-beam CT. The phantom has l marking points, where the distance between the jth marking point and a selected marking point in the spatial coordinate system is D real ; Among them, the spatial coordinate system takes the center of the turntable as its origin, which is the same as the axis vector of the world coordinate system;

[0009] Step 2: Determine the actual projection data (v ij ,w ij ), where the detector local coordinate system takes the detector center as the origin, and its three coordinate axes rotate clockwise around the x, y, and z axes by γ, θ respectively. y ,θ z Spend;

[0010] Step 3: Obtain the predicted projection data (v i (ζ j ),w i (ζ j ));

[0011] Step 4: The actual projection data (v ij ,w ij ) and the predicted projection data (v i (ζ j ),w i (ζ j )) Set as the objective function fun, and use the optimization method to obtain the optimal geometric parameter vector of the jth marker point. The geometric parameter vector includes the distance D between the jth marker point and the selected marker point in the spatial coordinate system. rec , and the distance sod between the ray source and the center of the turntable;

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

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

[0014]

[0015] In the formula, ζ j represents the geometric parameter vector of the jth marker point, ζ up ,ζ down Respectively represent ζ 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 of step 4 specifically includes:

[0017] Step 41, transform 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 optimal geometric parameter vector best_sol, assign fun(x0) to fun(current_sol) and the historical optimal geometric parameter vector fun(best_sol), fun(current_sol) = current_cost, fun(best_sol) = best_cost;

[0019] Step 43, using current_sol in step 42 as the initial value vector, using the first simplex method to obtain new solutions current_sol and (current_cost), if fun(current_sol) is less than best_cost, assign current_sol to best_sol, and assign fun(current_sol) to best_cost;

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

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

[0022] Step 431, create n linearly independent vertices X according to the initial value vector p =X0+Q p *e p , Qp represents the step length corresponding to the pth vertex, e p represents the n-dimensional unit vector corresponding to the p-th vertex. If 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 parameters is limited to make the convergence faster.

[0023] Step 432, when the optimization method loops to the preset number of cycles of the first simplex method, the second simplex method is interspersed to find a solution that meets the requirements as the initial value vector of the next first simplex method; wherein the preset number of cycles of the second simplex method is less than that of the first simplex method.

[0024] Furthermore, the marking point is the center of mass corresponding to the preset ball.

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

[0026] The phantom is set at the center of the turntable of the actual cone-beam CT. The phantom has l marking points, where the distance between the jth marking point and a selected marking point in the spatial coordinate system is D real ; Among them, the spatial coordinate system takes the center of the turntable as its origin, which is the same as the axis vector of the world coordinate system;

[0027] The actual projection data acquisition unit is used to determine the actual projection data (v ij ,w ij ), where the detector local coordinate system takes the detector center as the origin, and its three coordinate axes rotate clockwise around the x, y, and z axes by γ, θ respectively. y ,θ z Spend;;

[0028] The predicted projection data acquisition unit is used to acquire the predicted projection data (v i (ζ j ),w i (ζ j ));

[0029] The geometric parameter optimization unit is used to obtain the actual projection data (v ij ,w ij ) and the predicted projection data (v i (ζ j ),w i (ζ j )) Set as the objective function fun, and use the optimization method to obtain the optimal geometric parameter vector of the jth marker point. The geometric parameter vector includes the distance D between the jth marker point and the selected marker point in the spatial coordinate system. rec , and the distance sod between the ray source and the center of the turntable;

[0030] The distance correction unit is used to correct the D real , and the D of the geometric parameter optimization unit rec and sod, obtain the corrected distance sod between the jth marker and the selected marker in the spatial coordinate system real .

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

[0032]

[0033] In the formula, ζ j represents the geometric parameter vector of the jth marker point, ζ up ,ζ down Respectively represent ζ 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 of the geometric parameter optimization unit specifically includes:

[0035] The initial value vector setting subunit 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 subunit is used to calculate the objective function fun(X0) at X0, assign X0 to the current geometric parameter vector current_sol and the historical optimal geometric parameter vector best_sol, assign fun(x0) to fun(current_sol) and the historical optimal geometric parameter vector fun(best_sol), fun(current_sol) = current_cost, fun(best_sol) = best_cost;

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

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

[0039] Furthermore, the simplex method subunit specifically includes:

[0040] Vertex creation module, which is used to create n linearly independent vertices X according to the initial value vector p =X0+Q p *e p , Q p represents the step length corresponding to the pth vertex, e p represents the n-dimensional unit vector corresponding to the p-th vertex. If 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 parameters is limited to make the convergence faster.

[0041] A loop number judgment module is used to intersperse the use of the second simplex method when the optimization method loops to the preset loop number of the first simplex method 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.

[0042] Furthermore, the marking point is the center of mass corresponding to the preset ball.

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

[0044] 1. The present invention uses the Nelder–Mead method repeatedly, that is, after using the Nelder–Mead method continuously, a new solution is randomly found in its neighborhood as the initial point of the next Nelder–Mead method. Therefore, the present invention can overcome the problem of falling into the local minimum due to a bad initial value vector.

[0045] 2. The present invention uses the large and small simplex methods to search alternately and uses a penalty function to keep the simplex within a reasonable range. Therefore, the present invention can solve the problem of often taking too long to find a solution that meets the expectations and its stability.

[0046] 3. After searching for a relatively good solution, the present invention limits the search range of certain parameters that have a great influence on the objective function, which can make the convergence faster and limit the parameter changes when creating the vertices of the simplex. Therefore, the present invention can obtain a better solution and converge faster. BRIEF DESCRIPTION OF THE DRAWINGS

[0047] Figure 1 Schematic diagram of a phantom according to an embodiment of the present invention. DETAILED DESCRIPTION

[0048] In the drawings, the same or similar reference numerals are used to represent the same or similar elements or elements with the same or similar functions. The embodiments of the present invention are described in detail below in conjunction with the drawings.

[0049] In the description of the present invention, the terms "center", "longitudinal", "lateral", "front", "back", "left", "right", "vertical", "horizontal", "top", "bottom", "inside", "outside" and the like indicate directions or positional relationships based on the directions or positional relationships shown in the accompanying drawings, and are only for the convenience of describing the present invention and simplifying the description, rather than indicating or implying that the device or element referred to must have a specific direction, be constructed and operated in a specific direction, and therefore should not be understood as limiting the scope of protection of the present invention.

[0050] The method for correcting geometric parameters of a high-precision industrial CT imaging system provided by an embodiment of the present invention includes:

[0051] Step 1: Set the phantom at the center of the turntable of the actual cone-beam CT. The phantom has l marking points, where the distance between the jth marking point and a selected marking point in the spatial coordinate system is D real .like Figure 1 As shown, this embodiment uses a model composed of two forest balls (rubies) with a known distance, that is, l = 2. The marking point is the center of mass corresponding to the pre-set small ball, that is, D real is a known parameter and is the distance between the centers of mass of the two forest balls.

[0052] The spatial coordinate system takes the center of the turntable as its origin, which is the same as the axis vector of the world coordinate system. The axis vectors of the world coordinate system include the x, y, and z axes. The x axis is the axis of the line connecting the ray source and the detector, the z axis is the axis perpendicular to the horizontal plane, and the y axis is perpendicular to the x and z axes. Rotating around the x axis yields γ, and rotating around the y axis yields θ. y , rotate around the z axis to obtain θ z .

[0053] Step 2: Determine the actual projection data (v ij ,w ij ). The detector local coordinate system takes the detector center as the origin, and its three coordinate axes rotate clockwise around the x, y, and z axes by γ, θ respectively. y ,θ z Spend.

[0054] Step 3: Obtain the predicted projection data (v i (ζ j ),w i (ζ j )).

[0055] Step 4: The actual projection data (v ij ,w ij ) and the predicted projection data (v i (ζ j ),w i (ζ j )) Set as the objective function fun, and use the optimization method to obtain the optimal geometric parameter vector of the jth marker point. The geometric parameter vector includes the distance D between the jth marker point and the selected marker point in the spatial coordinate system. rec , the distance between the ray source and the center of the turntable sod, the distance between the straight line through the center of the turntable and the intersection of the detection and the ray source sdd, the deviation of the detector center from the ideal position coordinate in the detector plane uv (offset u ,offset v ), the rotation angle of the detector relative to the world coordinate system (θ y ,θ z ,γ), and the coordinates of the jth marker point in the spatial coordinate system (x j ,y j ,z j ).

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

[0057] In one embodiment, the optimization method takes minimizing fun as the optimization goal, as shown in formula (1):

[0058]

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

[0060] Of course, the optimization method can also achieve other optimization goals besides minimization according to actual needs.

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

[0062] Step 41, transform the n*1 dimensional geometric parameter vector ζ j The initial value vector of is set to X0. Where n represents ζ j The number of parameters included, such as ζ j When it 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 optimal geometric parameter vector best_sol, assign fun(x0) to fun(current_sol) and the historical optimal geometric parameter vector fun(best_sol), fun(current_sol) = current_cost, fun(best_sol) = best_cost.

[0064] Step 43, use current_sol in step 42 as the initial value vector, and use the first simplex method to obtain new solutions current_sol and fun(current_cost). If fun(current_sol) is less than best_cost, assign current_sol to best_sol and assign fun(current_sol) to best_cost.

[0065] Step 44, use the improved neighborhood search method at current_sol to search for a random point in a neighborhood as a new initial value vector new_sol: if fun(new_sol) is less than current_cost, assign new_sol to current_sol; if fun(new_sol) is not less than current_cost, and the current state of the loop has met the preset loop requirement, such as executing this operation every n times, assign new_sol to current_sol; if fun(new_sol) is not less than current_cost, and the current state of the loop has not met the preset loop requirement, 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 stronger ability to jump out of the local minimum solution before a better solution is found. When a better solution is found (that is, when best_sol is relatively small), certain key geometric parameters can be limited to randomly search the neighborhood range so that the objective function value of the updated solution is smaller, thereby obtaining a solution with a smaller objective function after using the simplex method.

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

[0068] Step 431, create n linearly independent vertices X according to the initial value vector p =X0+Q p *e p , Q p is a number randomly selected from a certain range, indicating the step length corresponding to the pth vertex, e pRepresents the n-dimensional unit vector corresponding to the p-th vertex, and the subsequent steps are consistent with the standard simplex method steps. If current_sol is less than the second preset threshold that can be obtained through multiple simplex method iterations, such as a pixel size, then the step size of key parameters (key parameters such as offset_x and γ) is limited to make the convergence faster, which can replace the neighborhood random search to a certain extent to jump out of the local minimum solution, and this method can make the previous simplex get a better point as the initial value vector, which is more efficient.

[0069] Step 432, when optimization method loops to the preset cycle number of the first simplex method, then intersperse and use the second simplex method, for example can be set to use the second simplex method every n times, the purpose of the second simplex method cycle setting here is to replace the neighborhood random search to a certain extent and make the search that jumps out of the local minimum solution, and can find the smaller solution of objective function value as the initial value vector of the first simplex method next time.Wherein, the preset cycle number of the second simplex method is less than the first simplex method.The present embodiment can be set with larger random step length to jump out of the local minimum solution at the second simplex method place, can set shorter step length before some parameters that affect the objective function larger as shown in the above-mentioned steps 431 when the initial value vector meets the requirements of the first simplex method, the second simplex method can take into account speed and can obtain a better solution relative to randomly carrying out random neighborhood.

[0070] In this embodiment, simplexes of random sizes are used continuously, which to a certain extent takes into account the ability to quickly find a relatively good solution and jump out of the local optimal solution, and combines a flexible neighborhood random search algorithm. When current_sol is less than one pixel, the solution is considered good enough. At this time, some parameters that have a greater impact on the objective function value can be limited.

[0071] The embodiment of the present invention further provides a high-precision correction system for geometric parameters of an industrial CT imaging system, which includes 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 set at the center of the turntable of the actual cone-beam CT. The phantom has l marking points, where the distance between the jth marking point and a selected marking point in the spatial coordinate system is D real ; Among them, the spatial coordinate system takes the center of the turntable as its origin, which is the same as the axis vector of the world coordinate system.

[0073] The actual projection data acquisition unit is used to determine the actual projection data (v ij ,w ij), where the detector local coordinate system takes the detector center as the origin, and its three coordinate axes rotate clockwise around the x, y, and z axes by γ, θ respectively. y ,θ z Spend.

[0074] The predicted projection data acquisition unit is used to acquire the predicted projection data (v i (ζ j ),w i (ζ j )).

[0075] The geometric parameter optimization unit is used to obtain the actual projection data (v ij ,w ij ) and the predicted projection data (v i (ζ j ),w i (ζ j )) Set as the objective function fun, and use the optimization method to obtain the optimal geometric parameter vector of the jth marker point. The geometric parameter vector includes the distance D between the jth marker point and the selected marker point in the spatial coordinate system. rec , and the distance sod between the ray source and the center of the turntable.

[0076] The distance correction unit is used to correct the D real , and the D of the geometric parameter optimization unit rec and sod, obtain the corrected distance sod between the jth marker and the selected marker 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 subunit 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, assign X0 to the current geometric parameter vector current_sol and the historical optimal geometric parameter vector best_sol, and assign fun(x0) to fun(current_sol) and the historical optimal geometric parameter vector fun(best_sol), fun(current_sol) = current_cost, fun(best_sol) = best_cost.

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

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

[0082] In one embodiment, the simplex method subunit specifically includes a vertex creation module and a loop number determination module, wherein:

[0083] The vertex creation module is used to create n linearly independent vertices X according to the initial value vector. p =X0+Q p *e p , Q p represents the step length corresponding to the pth vertex. p can be set to a random number in a reasonable interval. p Represents the n-dimensional unit vector corresponding to the p-th vertex. If current_sol is less than a second preset threshold that can be obtained through multiple simplex method iterations, the step size of the key parameters is limited to make convergence faster.

[0084] The loop number judgment module is used to intersperse the use of the second simplex method when the optimization method loops to the preset loop number of the first simplex method 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.

[0085] The advantages of the present invention are: first, a solution that meets the requirements can be found in a very short time and sufficiently stable. Second, the dependency of the initial value vector can be minimized as much as possible after determining an acceptable penalty function range. Not only is the motif required by this method very simple and easy to obtain, but the required sampling angle can also restore a relatively accurate solution in as few cases as possible.

[0086] The feasibility of this method is verified by simulation experiments. In the simulated circular trajectory scanning, the detector resolution is 1536*1536, the detector unit 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 ​​balls with a radius of 2mm are 25mm apart horizontally and 20mm apart vertically. When projected, the number of voxels is 1536*1536*1351, and the image size is 46.1612*46.1612*40.6829. The units of the geometric parameters of the distance in the table below are mm, and the units of the geometric parameters of the angle are 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 projections of the 60 centroids and the predicted projections of the centroids, which can be regarded as the Euclidean distance between them. At this time, the objective function is 0.0560. At this time, the distance error is about 0.2 pixels, which can be considered an acceptable solution.

[0092]

[0093] In the case that the true value is consistent with the above, one thousand sets of initial value vectors are randomly selected from the normal distribution in the interval 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]. The threshold is selected at 0.0561 to end the simplex and the large loop. It takes a total of 1906 seconds, and the average time to find a solution is 1.9 seconds each time.

[0094] The present invention proposes an optimization algorithm designed for the objective function used in the method to further search for the global minimum solution, namely the improved simplex method. Because the simplex method continuously uses simplex methods with different iteration times and simplex sizes, and a more flexible penalty function setting, the ability and speed of searching for the global optimal solution are improved on the basis of the simplex method. The present invention intersperses a small amount of random searches with restricted neighborhoods depending on the situation during the iteration process, so that the minimum solution can be more effectively jumped out while receiving the searched solution as good as possible, so that it will have a better effect as the initial value vector of the next simplex method.

[0095] Finally, it should be pointed out that the above embodiments are only used to illustrate the technical solutions of the present invention, rather than to limit them. Those skilled in the art should understand that the technical solutions described in the above embodiments may be modified, or some of the technical features thereof may be replaced by equivalents; these modifications or replacements do not deviate the essence of the corresponding technical solutions from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. A high-precision method for correcting geometric parameters of an industrial CT imaging system, characterized in that: include: Step 1: Set the phantom at the center of the turntable of the actual cone-beam CT. The phantom has l marking points, where the distance between the jth marking point and a selected marking point in the spatial coordinate system is D real ; Among them, the spatial coordinate system takes the center of the turntable as its origin, which is the same as the axis vector of the world coordinate system; Step 2: Determine the actual projection data (v ij ,w ij ), where the detector local coordinate system takes the detector center as the origin, and its three coordinate axes rotate clockwise around the x, y, and z axes by γ, θ respectively. y ,θ z Spend; Step 3: Obtain the predicted projection data (v i (ζ j ),w i (ζ j )); Step 4: The actual projection data (v ij ,w ij ) and the predicted projection data (v i (ζ j ),w i (ζ j )) Set as the objective function fun, and use the optimization method to obtain the optimal geometric parameter vector of the jth marker point. The geometric parameter vector includes the distance D between the jth marker point and the selected marker point in the spatial coordinate system. rec , and the distance sod between the ray source and the center of the turntable; Step 5, use D from step 1 real , and D in step 4 rec and sod, obtain the corrected distance sod between the jth marker and the selected marker in the spatial coordinate system real .

2. The method for correcting geometric parameters of a high-precision industrial CT imaging system according to claim 1, characterized in that: The optimization method takes minimizing fun as the optimization goal, as shown in formula (1): In the formula, ζ j represents the geometric parameter vector of the jth marker point, ζ up ,ζ down Respectively represent ζ j The upper and lower bounds of F(ζ up ,ζ down ) represents the penalty function term, and K represents the number of sampling angles.

3. The method for correcting geometric parameters of a high-precision industrial CT imaging system according to claim 1 or 2, characterized in that: The optimization method of step 4 specifically includes: Step 41, transform 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; Step 42, calculate the objective function fun(X0) at X0, assign X0 to the current geometric parameter vector current_sol and the historical optimal geometric parameter vector best_sol, assign fun(x0) to fun(current_sol) and the historical optimal geometric parameter vector fun(best_sol), fun(current_sol) = current_cost, fun(best_sol) = best_cost; Step 43, using current_sol in step 42 as the initial value vector, using the first simplex method to obtain new solutions current_sol and fun(current_cost); if fun(current_sol) is less than best_cost, assign current_sol to best_sol, and assign fun(current_sol) to best_cost; Step 44, at current_sol, use the improved neighborhood search method to search for a random point in a neighborhood as a new initial value vector new_sol: if fun(new_sol) is less than current_cost, assign new_sol to current_sol; if fun(new_sol) is not less than current_cost, and the current state of the loop has met the preset loop requirement, assign new_sol to current_sol; if fun(new_sol is not less than current_cost, and the current state of the loop has not met the preset loop requirement, execute step 43 until best_sol is less than the first preset threshold.

4. The method for correcting geometric parameters of a high-precision industrial CT imaging system according to claim 3, characterized in that: The first simplex method of step 43 specifically includes: Step 431, create n linearly independent vertices X according to the initial value vector p =X0+Q p *e p , Q p represents the step length corresponding to the pth vertex, e p represents the n-dimensional unit vector corresponding to the p-th vertex. If 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 parameters is limited to make the convergence faster. Step 432, when the optimization method loops to the preset number of cycles of the first simplex method, the second simplex method is interspersed to find a solution that meets the requirements as the initial value vector of the next first simplex method; wherein the preset number of cycles of the second simplex method is less than that of the first simplex method.

5. The method for calibrating geometric parameters of a high-precision industrial CT imaging system according to claim 1, characterized in that: The marking point is the center of mass of the pre-set ball.

6. A high-precision correction system for geometric parameters of industrial CT imaging systems, characterized in that: include: The phantom is set at the center of the turntable of the actual cone-beam CT. The phantom has l marking points, where the distance between the jth marking point and a selected marking point in the spatial coordinate system is D real ; Among them, the spatial coordinate system takes the center of the turntable as its origin, which is the same as the axis vector of the world coordinate system; The actual projection data acquisition unit is used to determine the actual projection data (v ij ,w ij ), where the detector local coordinate system takes the detector center as the origin, and its three coordinate axes rotate clockwise around the x, y, and z axes by γ, θ respectively. y ,θ z Spend;; The predicted projection data acquisition unit is used to acquire the predicted projection data (v i (ζ j ),w i (ζ j )); The geometric parameter optimization unit is used to obtain the actual projection data (v ij ,w ij ) and the predicted projection data (v i (ζ j ),w i (ζ j )) Set as the objective function fun, and use the optimization method to obtain the optimal geometric parameter vector of the jth marker point. The geometric parameter vector includes the distance D between the jth marker point and the selected marker point in the spatial coordinate system. rec , and the distance sod between the ray source and the center of the turntable; The distance correction unit is used to correct the D real , and the D of the geometric parameter optimization unit rec and sod, obtain the corrected distance sod between the jth marker and the selected marker in the spatial coordinate system real .

7. The high-precision industrial CT imaging system geometric parameter correction system according to claim 6, characterized in that: The optimization method takes minimizing fun as the optimization goal, as shown in formula (1): In the formula, ζ j represents the geometric parameter vector of the jth marker point, ζ up ,ζ down Respectively represent ζ j The upper and lower bounds of F(ζ up ,ζ down ) represents the penalty function term, and K represents the number of sampling angles.

8. The high-precision industrial CT imaging system geometric parameter correction system according to claim 6 or 7, characterized in that: The optimization method of the geometric parameter optimization unit specifically includes: The initial value vector setting subunit 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; The initial value calculation subunit is used to calculate the objective function fun(X0) at X0, assign X0 to the current geometric parameter vector current_sol and the historical optimal geometric parameter vector best_sol, assign fun(x0 to fun(current_sol) and the historical optimal geometric parameter vector fun(best_sol), fun(current_sol) = current_cost, fun(best_sol) = best_cost; The simplex method subunit is used to use the current_sol of the initial value calculation subunit as the initial value vector, and use the first simplex method to obtain new solutions current_sol and fun(current_cost). If fun(current_sol is less than best_cost, current_sol is assigned to best_sol, and fun(current_sol) is assigned to best_cost. The neighborhood search subunit is used to use the improved neighborhood search method at the current_sol to search for a random point in a neighborhood as a new initial value vector new_sol: if fun(new_sol) is less than current_cost, new_sol is assigned to current_sol; if fun(new_sol) is not less than current_cost, and the current state of the loop has met the preset loop requirements, new_sol is assigned to current_sol; if fun(new_sol) is not less than current_cost, and the current state of the loop has not met the preset loop requirements, the best_sol obtained by the simplex method subunit is less than the first preset threshold.

9. The high-precision industrial CT imaging system geometric parameter correction system according to claim 8, characterized in that: The simplex method subunits specifically include: Vertex creation module, which is used to create n linearly independent vertices X according to the initial value vector p =X0+Q p *e p , Q p represents the step length corresponding to the pth vertex, e p represents the n-dimensional unit vector corresponding to the p-th vertex. If 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 parameters is limited to make the convergence faster. A loop number judgment module is used to intersperse the use of the second simplex method when the optimization method loops to the preset loop number of the first simplex method 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.

10. The high-precision industrial CT imaging system geometric parameter correction system according to claim 6, characterized in that: The marking point is the center of mass of the pre-set ball.

Images