CN104764756B - The scaling method of cone-beam CT imaging system - Google Patents

Abstract

The invention discloses a calibration method for a cone-beam CT imaging system. The projected coordinates of 2N steel ball centers, and respectively fit the trajectory T ₁ and trajectory T ₂ of the ellipse where the projected coordinates of the N steel ball centers of each group are located; obtain the origin O and intersection point D of the i system and obtain the angles η, A ₁ , Four coordinates of B ₁ , A ₂ , and B ₂ ; construct the first objective function f(θ*, w _y *, w _z *) to calculate the angle θ and in the i system, calculate the origin O _w coordinate of the W coordinate system ( w _y , w _z ) and light source P _s coordinates (s _y , s _z ); construct the second objective function and calculate the angle Construct the third objective function H(t) and find the angle t; calculate the coordinates of the light source P _s in the W system and origin O coordinates ( ). The invention can calibrate the CBCT system of arbitrary complex orbits without knowing any geometric information of the system in advance, and can obtain all the geometric parameters of the projected position by using one projection, and has strong universality and pertinence.

Description

Calibration Method of Cone Beam CT Imaging System

technical field

The invention relates to the technical field of X-ray computed tomography imaging, in particular to a calibration method for a cone-beam CT imaging system.

Background technique

Cone beam CT (Cone Beam Computed Tomography, CBCT) technology is to use the cone beam X light source to rotate around the three-dimensional object at various angles, obtain the projection of the object at multiple angles, and then use the reconstruction algorithm of the cone beam projection map to reconstruct the three-dimensional object Image technology. Its principle is mainly to use the different characteristics of the attenuation coefficients of different components of the object on X-rays, and obtain a projection image that can reflect the characteristics of the object components on the detector behind the object.

The geometric parameters of the CBCT system refer to the geometric parameters that describe the angle and position between the X light source, object and detector of the system. These geometric parameters will affect the projected image of the object and will be used as important information in the reconstruction algorithm. Therefore, the geometric parameter calibration of CBCT system is of great significance to improve the reconstructed image quality and prevent image artifacts.

In the process of practical application, the scanning trajectory of the X light source in the CBCT system is not necessarily a complete circle. For example, in the related art, the light source of the C-arm CT system is fixed on the C-shaped rotating arm, and the doctor can operate on the patient at the same time. , observe the imaging situation in time. The scanning trajectory of this C-arm CT system is not a complete circle. In addition, due to the effect of gravity, the center of the trajectory will deviate from the center of the circle, causing the distance from the X-ray source to the detector to change. This scanning along an irregular trajectory The CBCT system also needs to accurately calibrate its geometric parameters.

Contents of the invention

The present invention aims to solve one of the technical problems in the related art at least to a certain extent. Therefore, an object of the present invention is to propose a calibration method for a cone-beam CT imaging system, which can calibrate a CBCT system with any complex orbit, and has strong universality.

In order to achieve the above object, the calibration method of the cone beam CT imaging system in the embodiment of the first aspect of the present invention includes the following steps: constructing a phantom, wherein the phantom is a calibration reference object, and the phantom is generally a plastic cylinder Two sets of steel balls are inlaid on the upper and lower bottom surfaces of the mold body, each group includes N (N is not less than 5) steel balls, the diameter of the bottom surface of the mold body is recorded as 2r, and the center of the two groups of steel balls determines The distance between two planes is recorded as 2l; The values of all pixels in the system are weighted to calculate the weighted average coordinate value of the projection of each steel ball to obtain the projected coordinates of the 2N steel ball centers, and respectively fit the trajectory of the ellipse where the projected coordinates of the N steel ball centers of each group are located T ₁ and track T ₂ ; the projected coordinates of the N steel ball centers of each group are respectively ordered on the corresponding track, and obtained The origin O and the intersection point D and the coordinates of the calculated angle η and A ₁ , B ₁ , A ₂ , and B ₂ under the system, wherein, the projected coordinates of the 2N steel ball centers are all in the Expressed below, the A ₁ and B ₁ are the projections of P ₁ and Q ₁ on the OD respectively, and the A ₂ and B ₂ are the projections of P ₂ and Q ₂ on the OD respectively, which is called the rectangle P ₁ P ₂ Q ₁ Q ₂ (length 2l, width 2r) is the detector around the under the After rotating the angle θ, the phantom is under the The section cut by the plane, the angle η is the detector around the rotate the angle θ and around the under the Rotation angle after, around the under the angle of rotation; in the In the plane, according to the angle η, the origin O, the intersection point D, Q ₁ coordinates (s ₁ , t ₁ ), P ₂ coordinates (s ₂ , t ₂ ), A ₁ coordinates A ₂ coordinates B ₁ coordinates B ₂ coordinates Construct the first objective function f(θ*,w _y *,w _z *), and use it to calculate the angle θ, in the The origin O _w coordinates (w _y , w _z ) and the light source P _s coordinates (s _y , s _z ) in the plane, wherein, the origin O _w is also the center of the phantom; according to the The projection of any point on the two bottom circles of the phantom should be respectively located _on the geometric relationship between the trajectory T1 and the trajectory T2 to construct the _second objective function and calculate the angle According to the origin O, the angle θ, the angle η, the angle and in the The positional relationship between the origin O _w coordinates (w _y , w _z ) and the light source P _s coordinates (s _y , s _z ) in the plane and the steel ball projection center coordinates constructs a third objective function H(t) and obtains out angle t; and according to the origin O, the angle θ, the angle η, the angle The angle t and the The origin O _w coordinates (w _y , w _z ) and the light source P _s coordinates (s _y , s _z ) in the plane are calculated at The coordinates of the light source P _s in the system and origin O coordinates

The calibration method of the cone beam CT imaging system according to the embodiment of the present invention has the following beneficial effects: (1) It can calibrate the CBCT system of any complex orbit, and has strong universality. (2) All the geometric parameters of the projection position can be obtained by using one projection, so it is very targeted for the situation that the scanning angle of the CT system is incomplete. (3) It does not need to know the geometric information of any system in advance, and can be used for geometric systems with poor mechanical precision.

Further, in one embodiment of the present invention, the The system is a world coordinate system that is stationary relative to the surrounding environment, and the With the center of the phantom as the origin O _w , the The axis of the phantom is Before starting projection, the tied perpendicular to the The ray is said under the according to the said and the right-hand rule; the is the virtual detector coordinate system describing the ideal state of the detector, the The point of intersection of the connection line between the light source P _s and the origin O _w and the detector plane is the origin O, the and said same direction, define P _s O with the Determine the plane for the plane, the through the origin O and perpendicular to the according to the said and the right-hand rule determine; and the system is the coordinate system based on the detector itself, the With the origin O as the origin, the and said parallel to the row and column lines of the detectors, the perpendicular to the detector plane, where, defining the said said respectively with the said said The corresponding coincident state is the initial state of the detector.

Further, in one embodiment of the present invention, from the initial state of the detector to any state of the detector includes the following steps: the detector moves around the rotation, with the form the angle θ, where the value of the angle θ is in accordance with the is defined in a right-hand rule manner for the positive direction; the detector orbits the rotation, with the form the angle Here, the angle value as described in is defined in a right-hand rule manner for the positive direction; and the detector orbits the rotation, with the Form the angle η, where the value of the angle η is according to the Defined in the right-hand rule manner for the positive direction.

Further, in one embodiment of the present invention, the trajectory T ₁ and trajectory T ₂ of the ellipse where the projected coordinates of the centers of the N steel balls of each group are respectively fitted according to the following equations:

Among them, a ₁ , b ₁ , c ₁ , and a ₂ , b ₂ , c ₂ , Projected by the centers of the N steel balls in each group on the The coordinates of the system are fitted.

Further, in one embodiment of the present invention, the projected coordinates of the center projection coordinates of the N steel balls in each group are respectively ordered on the corresponding trajectory, and obtained The origin O and the intersection point D under the system are specifically as follows: the projected coordinates of the N steel ball centers of each group are respectively ordered on the corresponding trajectory; the origin O is obtained by connecting E ₁ F ₂ and E ₂ F ₁ , wherein, The E ₁ and the E ₂ are respectively the steel ball center projection with the largest ordinate and the steel ball center projection with the smallest ordinate on the trajectory T ₁ , and the F ₁ and the F ₂ are respectively the projections of the steel ball center on the trajectory T ₂ . The steel ball center projection with the largest ordinate and the steel ball center projection with the smallest ordinate; and connecting E ₁ F ₁ and E ₂ F ₂ , the intersection of the extension lines of E ₁ F ₁ and E ₂ F ₂ in space as the all The intersection point D, wherein, the OD and the The included angle is the angle η.

Further, in one embodiment of the present invention, the angle θ is calculated according to the following steps, in the The origin O _w coordinates (w _y , w _z ) and the light source P _s coordinates (s _y , s _z ) in the plane: in the In the plane, for the undetermined origin O _w * coordinates (w _y *, w _z *), according to the length and width of the rectangle P ₁ P ₂ Q ₁ Q ₂ to obtain the corresponding undetermined rectangle vertices P ₁ *, P ₂ * , Q ₁ *, Q ₂ *, connect the origin O and the origin O _w * to obtain a straight line OO _w *; undetermined angle θ*, pass the intersection D to make the The parallel line DP _s * and said straight line OO _w * intersect at point P _s *, connecting A ₁ P ₁ *, B ₁ Q ₁ *, A ₂ P ₂ *, B ₂ Q ₂ * with said straight line OO _w * intersects at Construct the first objective function f(θ*, w _y *, w _z *) to represent the P _s *, degree of dispersion, wherein the greater the first objective function f(θ*, w _y *, w _z *), the greater the P _s *, The greater the degree of dispersion, when the first objective function f(θ*,w _y *,w _z *)=0, it means that the P _s *, five-point coincidence; calculation The optimization results to obtain the angle θ and the The origin O _w coordinates (w _y , w _z ) in the plane, namely: θ=θ ₀ ; take θ ₀ , Corresponding to the obtained P _s *, The average coordinates of the five points are the light source P _s coordinates (s _y , s _z ).

Further, in one embodiment of the present invention, or Wherein, λ is generally 1 or 2.

Further, in one embodiment of the present invention, the angle is calculated according to the following steps Angle to be determined According to the obtained angle η, the angle θ, in the The origin O _w coordinates (w _y , w _z ), the light source P _s coordinates (s _y , s _z ) and the undetermined angle in the plane Compute the two base circles of the phantom in the Trajectories T ₁ *, T ₂ * of the system; on the trajectories T ₁ *, T ₂ *, respectively take points (u ₁ , v ₁ ) and (u ₂ , v ₂ ), the constructor To measure the degree of deviation of the points (u ₁ , v ₁ ) and (u ₂ , v ₂ ) from the track T ₁ and the track T ₂ respectively, wherein the The larger , the greater the degree of deviation between the points (u ₁ , v ₁ ) and (u ₂ , v ₂ ) and the track T ₁ and track T ₂ respectively; construct the second objective function to calculate the angle Wherein, M means that a total of M groups of points (u ₁ , v ₁ ) and (u ₂ , v ₂ ) on the trajectories T ₁ * and T ₂ * are taken,

Further, in one embodiment of the present invention,

Among them, a ₁ , b ₁ , c ₁ , and a ₂ , b ₂ , c ₂ , Projected by the centers of the N steel balls in each group on the The coordinates of the system are fitted.

Further, in one embodiment of the present invention, the angle t is calculated according to the following formula:

t=argmin{H(t*)}

in, The angle parameter t* describes the position of the steel ball on the bottom circle of the phantom, and the functions h ₁ (t*) and h ₂ (t*) are respectively when the angle η and the angle θ are known. , the origin O _w coordinates (w _y , w _z ), the light source P _s coordinates (s _y , s _z ) and the angle Under the premise of , the projection of the steel ball at the position corresponding to t* on the two bottom surfaces of the phantom and the Tie down the distance to the nearest steel ball projection center.

Further, in one embodiment of the present invention,

Description of drawings

1 is a flowchart of a calibration method of a cone beam CT imaging system according to an embodiment of the present invention;

2 is a schematic structural view of a phantom in a calibration method of a cone beam CT imaging system according to an embodiment of the present invention;

Fig. 3 is a schematic diagram of defining a coordinate system of a calibration method of a cone beam CT imaging system according to an embodiment of the present invention;

Fig. 4 is the calibration method of the cone beam CT imaging system according to an embodiment of the present invention, in which the detector rotates from the initial state of the detector Schematic diagram of the rotation of the rotation angle θ;

Fig. 5 is that the middle detector of the cone beam CT imaging system according to one embodiment of the present invention continues to circle Rotation angle Schematic diagram of the rotation;

Fig. 6 shows that the detector continues to circle around in the calibration method of the cone beam CT imaging system according to an embodiment of the present invention. Rotation schematic diagram of rotation angle η;

Fig. 7 shows that the detector continues to circle around in the calibration method of the cone beam CT imaging system according to an embodiment of the present invention. After rotation angle η Department and The schematic diagram of the coordinate system of the system;

Fig. 8 is a calibration method of a cone beam CT imaging system according to an embodiment of the present invention schematic diagram of the plane;

9 is a schematic diagram of the positioning of the origin O and the intersection D of the calibration method of the cone beam CT imaging system according to an embodiment of the present invention; and

Fig. 10 is a flowchart of a calibration method of a cone beam CT imaging system according to a specific embodiment of the present invention.

detailed description

Embodiments of the present invention are described in detail below, examples of which are shown in the drawings, wherein the same or similar reference numerals designate the same or similar elements or elements having the same or similar functions throughout. The embodiments described below by referring to the figures are exemplary and are intended to explain the present invention and should not be construed as limiting the present invention.

The following describes the calibration method of the cone beam CT imaging system according to the embodiment of the present invention with reference to the accompanying drawings.

Fig. 1 is a flowchart of a calibration method of a cone beam CT imaging system according to an embodiment of the present invention. As shown in Figure 1, the calibration method of the cone beam CT imaging system of the embodiment of the present invention includes the following steps:

S1, constructing the phantom, wherein, the phantom is a calibration reference object, the phantom is generally a plastic cylinder, two groups of steel balls are inlaid on the upper and lower bottom surfaces of the phantom, each group includes N (N is not less than 5) steel balls, the phantom The diameter of the body bottom circle is recorded as 2r, and the distance between the two planes determined by the centers of the two groups of steel balls is recorded as 2l.

Specifically, in one embodiment of the present invention, as shown in Figure 2, the centers of the two groups of steel balls are respectively arranged on the two bottom circles of the mold body, and N (for example, 6) steel balls in each group are placed equiangularly and evenly, The plane distance between the two bottom circles determined by the centers of two groups of steel balls is 2l, the diameter of the bottom circle determined by the centers of N steel balls in each group is 2r, and the radius of N, l, r and the steel balls themselves must meet the requirements to identify the steel balls on the projection map Under the premise of the center position, make the projected area as large as possible.

Further, for any CBCT system to be calibrated, first define the following three coordinate systems Tie, Department and Tie. As shown in Figure 3, The system is a world coordinate system that is stationary relative to the surrounding environment, Take the center of the phantom as the origin O _w , The axis of the phantom is (y axis, The forward and reverse are arbitrary), before starting the projection, tied perpendicular to A ray of (x-axis), under the (z-axis) according to and the right-hand rule determine. is the virtual detector coordinate system describing the ideal state of the detector, The intersection of the line connecting the light source P _s and the origin O _w with the detector plane is the origin O, (y-axis) and In the same direction, define P _s O and The determined plane is flat, (z-axis) through the origin O and perpendicular to (x-axis) according to and the right-hand rule determine, where, The positive direction of the light source P _s in The projection on Orientation when the half axis is positive. The system is the coordinate system based on the detector itself, With the origin O as the origin, (x-axis) and (y-axis) parallel to the row and column lines of the detectors, (z-axis) is perpendicular to the detector plane. Among them, define respectively with The state corresponding to coincidence is the initial state of the detector.

Further, in one embodiment of the present invention, from the initial state of the detector to any state of the detector includes the following steps:

S8, the detector orbits from the initial state of the detector rotate, with form the angle θ, where the value of the angle θ is according to Defined in the right-hand rule manner for the positive direction.

A schematic diagram of the rotation of the detector in step S8 is shown in FIG. 4 .

S9, the detector orbits from the initial state of the detector rotate, with form an angle Here, the angle value according to Defined in the right-hand rule manner for the positive direction.

A schematic diagram of the rotation of the detector in step S9 is shown in FIG. 5 .

S10, the detector rotates from the initial state of the detector to rotate, with form an angle η, where the value of the angle η is according to Defined in the right-hand rule manner for the positive direction.

The schematic diagram of the rotation of the detector in step S10 is shown in Figure 6, at this time Department and The system is shown in Figure 7.

It should be noted that, from step S8 to step S10, each rotation of the detector will make The direction of the two axes changes, in order to make Figure 4 to Figure 6 more understandable, in Figure 4 to Figure 6 refers to the detector before this rotation That is, the rotation angle θ, the angle The angle η is surrounded by

S2, projected by steel balls on The values of all pixels in the system are weighted to calculate the weighted average coordinate value of the projection of each steel ball to obtain the projected coordinates of 2N steel ball centers, and respectively fit the trajectory T ₁ of the ellipse where the projected coordinates of the N steel ball centers of each group are located and trajectory T ₂ .

Since the position of each group of N steel balls in space is on the bottom circle, the projected coordinates of the centers of the N steel balls in each group should be on the projection of the bottom circle, that is, on an ellipse. In one embodiment of the present invention, the locus T ₁ and the locus T ₂ of the ellipse where the projected coordinates of the N steel ball centers of each group are located can be respectively fitted according to the following equations:

Among them, a ₁ , b ₁ , c ₁ , and a ₂ , b ₂ , c ₂ , Projected by the centers of N steel balls in each group on The coordinates of the system are fitted.

S3, the projected coordinates of the center of N steel balls in each group are respectively ordered on the corresponding trajectory, and obtained The origin O, the intersection point D and the calculated angle η and the four-point coordinates of A ₁ , B ₁ , A ₂ , and B ₂ under the system, among which, the center projection coordinates of 2N steel balls are all in Expressed below, A ₁ and B ₁ are the projections of P ₁ and Q ₁ on OD respectively, A ₂ and B ₂ are the projections of P ₂ and Q ₂ on OD respectively, called rectangle P ₁ P ₂ Q ₁ Q ₂ (Length is 2l, width is 2r) for the detector around under the After rotating the angle θ, the phantom is under the The section cut by the plane, the angle η is the Rotation angle θ and around under the Rotation angle after, around under the The angle of rotation.

It should be noted that, as shown in Figure 8, the detector under the After rotating the angle θ, coincides with the line where OD is located and with The included angle is the angle θ, and the detector revolves around Rotation angle θ and around Rotation angle And around under the After shaft rotation angle η, OD and The included angle is angle η. In Fig. 8, ||P ₁ P ₂ ||=2l, ||P ₁ Q ₁ ||=2r, the straight line P _s D is parallel to And intersect with the detector plane at D, A ₁ B ₁ , A ₂ B ₂ are respectively obtained by intersecting the trajectory T ₁ and T ₂ with OD.

Further, in one embodiment of the present invention, the projected coordinates of the centers of N steel balls in each group are respectively ordered on the corresponding trajectories, and obtained The origin O and intersection point D under the system can be specifically:

S31. Respectively determine the order of the projected coordinates of the centers of the N steel balls in each group on the corresponding trajectories.

S32, connecting E ₁ F ₂ and E ₂ F ₁ to obtain the origin O, wherein E ₁ and E ₂ are respectively the steel ball center projection with the largest ordinate and the steel ball center projection with the smallest ordinate on the track T ₁ , F ₁ and F ₂ are respectively the projection of the steel ball center with the largest ordinate and the projection of the steel ball center with the smallest ordinate _on the trajectory T2.

E ₁ , E ₂ , F ₁ and F ₂ are shown in FIG. 3 .

S33, connect E ₁ F ₁ and E ₂ F ₂ , the intersection of E ₁ F ₁ and E ₂ F ₂ extended lines in space is taken as intersection D, where OD and The included angle is angle η.

A schematic diagram of positioning the origin O and the intersection point D in step S32 and step S33 is shown in FIG. 9 .

S4, in In the plane, according to angle η, origin O, intersection point D, Q ₁ coordinate (s ₁ , t ₁ ), P ₂ coordinate (s ₂ , t ₂ ), A ₁ coordinate A ₂ coordinates B ₁ coordinates B ₂ coordinates Construct the first objective function f(θ*,w _y *,w _z *), and use it to calculate the angle θ, in The origin O _w coordinates (w _y , w _z ) and the light source P _s coordinates (s _y , s _z ) in the plane, where the origin O _w is also the center of the phantom.

Specifically, in one embodiment of the present invention, the angle θ can be calculated according to the following steps, where The origin O _w coordinates (w _y , w _z ) and the light source P _s coordinates (s _y , s _z ) in the plane:

S41, at In the plane, for the undetermined origin O _w * coordinates (w _y *, w _z *), get the corresponding undetermined rectangle vertices P ₁ *, P ₂ *, Q according to the length and width of the rectangle P ₁ P ₂ Q ₁ Q _{2 1} *, Q ₂ *, connect the origin O and the origin O _w * to get the straight line OO _w *.

S42, undetermined angle θ*, pass the intersection point D for The parallel line DP _s * and the straight line OO _w * intersect at the point P _s *, connecting A ₁ P ₁ *, B ₁ Q ₁ *, A ₂ P ₂ *, B ₂ Q ₂ * intersect the straight line OO _w * at

S43. Construct the first objective function f(θ*, w _y *, w _z *) to represent P _s *, The degree of dispersion of , where the greater the first objective function f(θ*, w _y *, w _z *), the greater the P _s *, The greater the degree of dispersion, when the first objective function f(θ*,w _y *,w _z *)=0, it means P _s *, Five coincides.

S44, computing The optimization results of to obtain the angle θ and the The origin O _w coordinates (w _y , w _z ) in the plane, namely: θ=θ ₀ .

S45, take θ ₀ , Corresponding to the obtained P _s *, The average coordinates of the five points are the light source P _s coordinates (s _y , s _z ).

Specifically, in one embodiment of the present invention, s _z =||OD||sinθ, (s ₁ ,t ₁ )=(w _y -l,w _z -r), (s ₂ ,t ₂ )=(w _y +l,w _z +r). Further, according to P _s Q ₁ B ₁ , P _s P ₂ A ₂ are collinear, we can know: make

Further, in one embodiment of the present invention, when , O is the midpoint of B ₁ A ₂ , B ₁ A ₂ is parallel to Q ₁ P ₂ , Among them, 2r is the diameter of the bottom circle, and 2l is the distance between two planes determined by the centers of two groups of steel balls.

Further, in another embodiment of the present invention, when hour, w _y =ξw _z , in, or etc., where λ is generally 1 or 2.

S5, construct the second objective function according to the geometric relationship that the projection of any point on the _two bottom circles of the phantom should be located _on the trajectory T1 and the trajectory T2 respectively and calculate the angle

Further, in one embodiment of the present invention, the angle can be calculated according to the following steps

S51, angle to be determined According to the obtained angle η, angle θ, in The origin O _w coordinates (w _y , w _z ) in the plane, the light source P _s coordinates (s _y , s _z ) and the undetermined angle Compute the two base circles of the phantom at The trajectories T ₁ *, T ₂ * of the system.

S52. Take points (u ₁ , v ₁ ) and (u ₂ , v ₂ ) respectively on the trajectories T ₁ * and T ₂ *, and construct the function To measure the degree of deviation of points (u ₁ , v ₁ ) and (u ₂ , v ₂ ) from track T ₁ and track T ₂ respectively, where, The larger , the greater the degree of deviation of points (u ₁ , v ₁ ) and (u ₂ , v ₂ ) from track T ₁ and track T ₂ respectively.

Specifically, in one embodiment of the present invention,

Among them, a ₁ , b ₁ , c ₁ , and a ₂ , b ₂ , c ₂ , Projected by the centers of N steel balls in each group on The coordinates of the system are fitted.

S53, constructing a second objective function to calculate the angle Among them, M means that a total of M groups of points (u ₁ , v ₁ ) and (u ₂ , v ₂ ) on the trajectory T ₁ * and T ₂ * have been taken,

Among them, M is an integer greater than or equal to 1, j is an integer smaller than M, and the second objective function used to measure the angle _The deviation between the projected coordinates of the center _of the lower steel ball and the trajectory T1 and T2.

It should be noted, plane around When rotating -η, The positions of all points on the plane are also around Rotated -η. At this time, at next, new The three basis vectors in the system can be expressed as:

Furthermore, in an embodiment of the present invention, for any point K _1,2 on the two bottom circles respectively, at next, Among them, α _1,2 ∈ [0,360°), define and The intersection point of the plane is R _1,2 , then in next, and Solutions have to: Find out according to λ exist The coordinates under the system, and then use can get R _1,2 in the new Coordinates under the system the coordinates Rotate the angle η around the origin O and then substitute into the equations _of the fitted trajectory T1 and trajectory T2 respectively. Therefore, the equations _{of trajectory} T1 and trajectory _T2 can be obtained according to any number _of points _on the trajectory T1 and trajectory T2 Each parameter in .

S6, according to origin O, angle θ, angle η, angle and in The position relationship between the origin O _w coordinates (w _y , w _z ) in the plane, the light source P _s coordinates (s _y , s _z ) and the center coordinates of the steel ball projection constructs the third objective function H(t) and obtains the angle t.

where the angle t is under the Positive semi-axis at under the projection on the plane and under the The included angle in the positive direction, and the sum of the positive direction of the angle t under the Satisfy the right-hand rule, where, in the first projection of the motif with coincide.

Since the first projection of the phantom corresponds to the with coincide, so the first projection of the phantom corresponds to t=0 of the system. Further, in an embodiment of the present invention, the angle t can be calculated according to the following formula:

in, The angle parameter t* describes the position of the steel ball on the bottom circle of the phantom, and the functions h ₁ (t*), h ₂ (t*) are respectively at the known angle η, angle θ, and origin O _w coordinates ( w _y , w _z ), light source P _s coordinates (s _y , s _z ) and angle Under the premise of , the projection of the steel ball at the position corresponding to t* on the two bottom surfaces of the phantom and the Tie down the distance to the nearest steel ball projection center. Specifically, in one embodiment of the present invention, for After rotating the angle η around the origin O, to the nearest steel ball center distance, for After rotating the angle η around the origin O, to the nearest steel ball center distance, for P _s K _1,2 with The intersection point of the plane, K ₁ , 2 is any point on the two base circles respectively, in next, i is an integer greater than or equal to 0 and less than or equal to N-1. H(t) with For the period, determining the true value of H(t) needs to be combined with the angle t calibrated in the previous projection.

S7, according to origin O, angle θ, angle η, angle angle t and at The origin O _w coordinates (w _y , w _z ) and light source P _s coordinates (s _y , s _z ) in the plane are calculated at The coordinates of the light source P _s in the system and origin O coordinates

Further, in one embodiment of the present invention, The coordinates of the light source P _s in the system and the coordinates of the origin O It can be calculated according to the following formula:

The calibration method of the cone beam CT imaging system according to the embodiment of the present invention will be described below according to specific embodiments.

As shown in Figure 10, in a specific embodiment of the present invention, the calibration method of the cone beam CT imaging system includes the following steps:

S101, constructing a phantom.

Wherein, the test phantom takes N=6, l=40mm, r=65mm, the diameter of the steel ball itself is 2.5mm, and the total number of steel balls is 12. The detector of the system to be calibrated is 1920×1536 pixels, and the side length of the detector pixel unit is 0.127mm. For any CBCT system to be calibrated, define the coordinate system Tie, Tie, system, the rotation angle θ of the detector, η, and O to vertical winding The rotation angle t.

S102, calculate the weighted average coordinate value of the projection of each steel ball with the value of all pixels of the steel ball projection as weights, obtain the center coordinates of the 12 steel balls, and fit the trajectory T ₁ , T of the ellipse where the center coordinates of the 12 steel balls are located ₂ .

S103, order the centers of the 6 steel balls in each group, and obtain the origin O. At this time, the center coordinates of all 12 steel balls can be Under the expression.

S104, calculate the angle η, and obtain A ₁ B ₁ , A ₂ B ₂ from the intersection of trajectories T ₁ , T ₂ and OD.

S105, at In the plane, calculate the angle θ, the origin O _w coordinates (w _y , w _z ) corresponding to the angle θ, and the light source P _s coordinates (s _y , s _z ).

The angle θ error is about 0.01°, but the accuracy of the origin O _w coordinates (w _y , w _z ) and the light source P _s coordinates (s _y , s _z ) is not ideal (especially when θ→0,OD→∞) .

S106, performing a first calibration correction on the angle θ, the origin O _w coordinates (w _y , w _z ) and the light source P _s coordinates (s _y , s _z ).

If the angle θ>2° (2° is determined by the experimental conditions), then the origin The coordinates of are obtained (w _y ,w _z ). On the contrary, it is processed as follows: approximate the angle θ to 0. Assuming that the middle point of P ₁ Q ₁ is M ₁ , and the middle point of P ₂ Q ₂ is M ₂ , then according to "the projection of the intersection of two straight lines is equal to the intersection of their projections", the relative steel ball center can be connected in ellipses T ₁ and T ₂ And get the projection O ₁ , O ₂ of M ₁ , M ₂ by taking the intersection point of the connecting line (O ₁ , O ₂ are on the straight line OD).

Let ||OP _s ||/||OO _w ||=p', namely (s _y , s _z )=p'(w _y , w _z ), then ||O ₁ O ₂ ||/||M ₁ M ₂ ||＝p'/(p'-1), get:

S107, perform further calibration correction on the origin O _w coordinates (w _y , w _z ) and the light source P _s coordinates (s _y , s _z ).

Step S107 can also be repeated multiple times, the specific method is:

Take it within the range of (w _y ±0.2, w _z ±0.2) (0.2mm is the maximum error after the first calibration correction under the experimental conditions) Find the coordinates (w _y '±l,w _z '±r) of P ₁ P ₂ Q ₁ Q ₂ at this time, and then connect A ₁ P ₁ , B ₁ Q ₁ , A ₂ P ₂ , B ₂ Q ₂ and intersects with The objective function is express request variance), take w _y , w _z corresponding to It is the required (s _y , s _z ).

S108, get new system, and define subfunctions Define the objective function to measure Calculate the deviation between the projection of the center of the sphere and the trajectory of T ₁ and T ₂ under the angle, and calculate the angle

Pick: at this time α ₁ , α ₂ =10j, j=0, 1, 2...35 are on ellipses T ₁ , T ₂ .

S109, defining the objective function H(t) under the experimental conditions to calculate the angle t.

Pick:

S110, find out in P _s coordinates of the light source under the system and origin O coordinates

S111, the calibration ends.

The beneficial effects of the present invention are as follows: (1) The CBCT system of any complex orbit can be calibrated, and has strong universality. (2) All the geometric parameters of the projected position can be obtained by using one projection, which is very targeted for the situation that the scanning angle of the CT system is not complete. (3) It is not necessary to know the geometric information of any system in advance, and it can be used in geometric systems with poor mechanical precision, which is simple and convenient.

In the description of this specification, descriptions referring to the terms "one embodiment", "some embodiments", "example", "specific examples", or "some examples" mean that specific features described in connection with the embodiment or example , structure, material or characteristic is included in at least one embodiment or example of the present invention. In this specification, the schematic representations of the above terms are not necessarily directed to the same embodiment or example. Furthermore, the described specific features, structures, materials or characteristics may be combined in any suitable manner in any one or more embodiments or examples. In addition, those skilled in the art can combine and combine different embodiments or examples and features of different embodiments or examples described in this specification without conflicting with each other.

In addition, the terms "first" and "second" are used for descriptive purposes only, and cannot be interpreted as indicating or implying relative importance or implicitly specifying the quantity of indicated technical features. Thus, the features defined as "first" and "second" may explicitly or implicitly include at least one of these features. In the description of the present invention, "plurality" means at least two, such as two, three, etc., unless otherwise specifically defined.

Any process or method descriptions in flowcharts or otherwise described herein may be understood to represent modules, segments or portions of code comprising one or more executable instructions for implementing specific logical functions or steps of the process , and the scope of preferred embodiments of the invention includes alternative implementations in which functions may be performed out of the order shown or discussed, including substantially concurrently or in reverse order depending on the functions involved, which shall It is understood by those skilled in the art to which the embodiments of the present invention pertain.

The logic and/or steps represented in the flowcharts or otherwise described herein, for example, can be considered as a sequenced listing of executable instructions for implementing logical functions, can be embodied in any computer-readable medium, For use with instruction execution systems, devices, or devices (such as computer-based systems, systems including processors, or other systems that can fetch instructions from instruction execution systems, devices, or devices and execute instructions), or in conjunction with these instruction execution systems, devices or equipment used. For the purposes of this specification, a "computer-readable medium" may be any device that can contain, store, communicate, propagate or transmit a program for use in or in conjunction with an instruction execution system, device or device. More specific examples (non-exhaustive list) of computer-readable media include the following: electrical connection with one or more wires (electronic device), portable computer disk case (magnetic device), random access memory (RAM), Read Only Memory (ROM), Erasable and Editable Read Only Memory (EPROM or Flash Memory), Fiber Optic Devices, and Portable Compact Disc Read Only Memory (CDROM). In addition, the computer-readable medium may even be paper or other suitable medium on which the program can be printed, since the program can be read, for example, by optically scanning the paper or other medium, followed by editing, interpretation or other suitable processing if necessary. The program is processed electronically and stored in computer memory.

It should be understood that various parts of the present invention can be realized by hardware, software, firmware or their combination. In the above described embodiments, various steps or methods may be implemented by software or firmware stored in memory and executed by a suitable instruction execution system. For example, if implemented in hardware, as in another embodiment, it can be implemented by any one or combination of the following techniques known in the art: Discrete logic circuits, ASICs with suitable combinational logic gates, Programmable Gate Arrays (PGAs), Field Programmable Gate Arrays (FPGAs), etc.

Those of ordinary skill in the art can understand that all or part of the steps carried by the methods of the above embodiments can be completed by instructing related hardware through a program, and the program can be stored in a computer-readable storage medium. During execution, one or a combination of the steps of the method embodiments is included.

In addition, each functional unit in each embodiment of the present invention may be integrated into one processing module, each unit may exist separately physically, or two or more units may be integrated into one module. The above-mentioned integrated modules can be implemented in the form of hardware or in the form of software function modules. If the integrated modules are realized in the form of software function modules and sold or used as independent products, they can also be stored in a computer-readable storage medium.

The storage medium mentioned above may be a read-only memory, a magnetic disk or an optical disk, and the like. Although the embodiments of the present invention have been shown and described above, it can be understood that the above embodiments are exemplary and should not be construed as limiting the present invention, those skilled in the art can make the above-mentioned The embodiments are subject to changes, modifications, substitutions and variations.

Claims (10)

1. A calibration method of a cone beam CT imaging system is characterized by comprising the following steps:

building a die body, wherein the die body is a calibration reference object, the die body is generally a plastic cylinder, two groups of steel balls are respectively embedded in the upper bottom surface and the lower bottom surface of the die body, each group comprises N steel balls, the diameter of the bottom surface circle of the die body is recorded as 2r, the distance between two planes determined by the centers of the two groups of steel balls is recorded as 2l, and the number N of the steel balls in each group is not less than 5;

projected on the steel ballTaking the values of all pixels of the system as weights, solving a weighted average coordinate value for the projection of each steel ball to obtain the projection coordinates of the centers of 2N steel balls, and respectively fitting the trajectory T of the ellipse where the projection coordinates of the centers of N steel balls of each group are located₁And track T₂；

Respectively setting the central projection coordinates of the N steel balls in each group in a sequence on corresponding tracks, and obtaining an original point O and an intersection point D under the i system and calculation angles η and A₁、B₁、A₂、B₂Four-point coordinates, wherein the 2N steel ball central projection coordinates are all in theIs expressed under the line A₁、B₁Are respectively P₁、Q₁Projection on OD, said A₂、B₂Are respectively P₂、Q₂Projection on the OD, called rectangle P₁P₂Q₁Q₂For the detector aroundUnder the systemAfter rotating the angle theta, the model is y under the system of iⁱ-zⁱA cross section taken in a plane, said P₁，Q₁，P₂，Q₂Respectively is the rectangle P₁P₂Q₁Q₂The angle η being the angle around which the detector is atBy rotation of said angle theta and about saidUnder the systemRotation angleThen, wind around theUnder the systemAngle of rotation, wherein the origin O_wCoordinate (w)_y,w_z) Is the origin O_wY and z coordinates in the i system, light source P_sCoordinate(s)_y,s_z) Is a light source P_sY and z coordinates in the i system;

at said yⁱ-zⁱIn-plane, according to the angle η, the origin O, the intersection point D, Q₁Coordinate(s)₁,t₁)、P₂Coordinate(s)₂,t₂)、A₁Coordinates of the objectA₂Coordinates of the objectB₁Coordinates of the objectB₂Coordinates of the objectConstructing a first objective function f (theta, w)_y*,w_zAnd calculating the angle theta, y in this wayⁱ-zⁱThe origin O in the plane_wCoordinate (w)_y,w_z) And the light source P_sCoordinate(s)_y,s_z) Wherein the origin O_wAt the same time isThe center of the mold body;

the projections of any point on the two bottom surface circles of the die body are respectively positioned on the track T₁And said track T₂To construct a second objective functionAnd calculating the angle

According to the origin O, the angle theta, the angle η and the angleAnd at said yⁱ-zⁱThe origin O in the plane_wCoordinate (w)_y,w_z) And the light source P_sCoordinate(s)_y,s_z) Constructing a third objective function H (t) according to the position relation among the steel ball projection center coordinates, and solving an angle t, wherein the angle t describes the position of the steel ball on the bottom surface circle of the die body; and

according to the origin O, the angle theta, the angle η and the angleThe angle t and the angle yⁱ-zⁱThe origin O in the plane_wCoordinate (w)_y,w_z) And the light source P_sCoordinate(s)_y,s_z) Is calculated atLight source P under system_sCoordinates of the objectAnd origin O coordinateWherein,

the above-mentionedIs a world coordinate system that is stationary relative to the surrounding environment, saidUsing the center of the mold body as the origin O_wSaidIs about the axis of the mold bodyBefore starting projection, theIs perpendicular to theIs a ray ofThe above-mentionedUnder the systemAccording to the aboveThe above-mentionedAnd right hand rule determination;

the i is a virtual detector coordinate system describing the ideal state of the detector, and the i is the lightSource P_sAnd the origin O_wThe intersection point of the connecting line of (a) and the detector plane is the origin O, yⁱAnd saidIn the same direction, define P_sO and the yⁱThe determined plane is the yⁱ-zⁱPlane of said zⁱPassing through the origin O and perpendicular to the yⁱ，xⁱAccording to the yⁱZ toⁱAnd right hand rule determination; and

the above-mentionedIs based on the coordinate system of the detector itself, saidUsing the origin O as the originAnd saidParallel to the line of rows and columns of the detector, thePerpendicular to the detector plane, wherein theThe above-mentionedThe above-mentionedRespectively with said xⁱThe yⁱZ toⁱCorresponding to the state at the time of coincidence asThe initial state of the detector.

2. A method of calibrating a cone-beam CT imaging system according to claim 1, wherein the step of calibrating from an initial state of said detector to any state of said detector comprises the steps of:

the detector is wound around the detector from the initial state of the detectorRotate withForm said angle theta, where said angle theta value is in accordance withDefining a right-hand rule mode for the positive direction;

the detector is wound around the detector from the initial state of the detectorRotate withForm said angleHere, the angleValue is as followsDefining a right-hand rule mode for the positive direction; and

the detector is wound around the detector from the initial state of the detectorRotate withForming the angle η, where the angle η value is in accordance withDefined for the right-hand rule of the positive direction.

3. The calibration method of the cone-beam CT imaging system as claimed in claim 1, wherein the trajectory T of the ellipse where the central projection coordinates of the N steel balls in each group are located is fitted according to the following equation₁And track T₂：

Wherein, a₁、b₁、c₁、And a₂、b₂、c₂、The center of the N steel balls of each group is projected on theCoordinate fitting of the system.

4. The calibration method of the cone-beam CT imaging system according to claim 1, wherein the step of sequentially determining the central projection coordinates of the N steel balls in each set on the corresponding trajectory and obtaining the origin O and the intersection D under the i system comprises:

respectively setting the central projection coordinates of the N steel balls in each group in a corresponding track;

connection E₁F₂And E₂F₁Obtaining the origin O, wherein E₁And said E₂Respectively being the track T₁The center projection of the steel ball with the maximum vertical coordinate and the center projection of the steel ball with the minimum vertical coordinate are arranged, F₁And said F₂Respectively being the track T₂Projecting the center of the steel ball with the largest vertical coordinate and the center of the steel ball with the smallest vertical coordinate; and

connection E₁F₁And E₂F₂Said E is₁F₁And said E₂F₂The intersection point of the extension lines in the space is taken as the intersection point D, wherein the OD and theIs the angle η.

5. The method of claim 1, wherein the angle θ is calculated at the y-position according to the following stepsⁱ-zⁱThe origin O in the plane_wCoordinate (w)_y,w_z) And the light source P_sCoordinate(s)_y,s_z)：

At said yⁱ-zⁱIn the plane, for a pending origin O_wCoordinate (w)_y*,w_zAccording to said rectangle P)₁P₂Q₁Q₂The length and width of the rectangular object to be determined are obtained to obtain the corresponding vertex P of the rectangular object to be determined₁*、P₂*、Q₁*、Q₂Connecting the origin O and the origin O_wObtain a straight line OO_w*；

Angle θ to be determined, said y being taken as saidⁱParallel lines DP of_sAnd the straight line OO_wIntersect at point P_sIs connected to A₁P₁*、B₁Q₁*、A₂P₂*、B₂Q₂Respectively connected with the straight line OO_wIntersect at

Constructing a first objective function f (theta, w)_y*,w_zTo represent said P_s*、Wherein the first objective function f (θ, w)_y*,w_zGreater, said P_s*、When the first objective function f (theta, w) is larger than the first threshold value_y*,w_z0 denotes said P_s*、Overlapping five points;

computingTo obtain the angle theta and the angle yⁱ-zⁱThe origin O in the plane_wCoordinate (w)_y,w_z) Namely: theta is equal to theta₀(ii) a And

take theta₀,Corresponding to the obtained P_s*、The average coordinate of five points is the light source P_sCoordinate(s)_y,s_z)。

6. The method for calibrating a cone-beam CT imaging system as recited in claim 5, orWherein λ is generally 1 or 2.

7. The method of claim 1, wherein the angle is calculated according to the following steps

Angle to be determinedBased on the obtained angle η, the angle theta, and the yⁱ-zⁱThe origin O in the plane_wCoordinate (w)_y,w_z) The light source P_sCoordinate(s)_y,s_z) And angle to be determinedCalculating the two bottom surface circles of the die bodyTrajectory T of the system₁*、T₂*；

At the track T₁*、T₂On, take points (u) respectively₁,v₁) And (u)₂,v₂) Constructing a function q (u)₁,v₁,u₂,v₂,) To measure the point (u)₁,v₁) And (u)₂,v₂) Respectively with said track T₁And said track T₂Wherein, the degree of deviation ofThe larger the point (u)₁,v₁) And (u)₂,v₂) Respectively with said track T₁And track T₂The greater the degree of deviation; and

constructing a second objective functionTo calculate said angleWherein M represents that M groups are taken out in the track T₁*、T₂Said point (u) on₁,v₁) And (u)₂,v₂)，

8. The method of calibrating a cone-beam CT imaging system as recited in claim 7,

。

9. a method of calibrating a cone-beam CT imaging system according to claim 1, wherein said angle t is calculated according to the following formula:

t＝argmin{H(t*)}

wherein,t is the angle to be determined of said angle t, function h₁(t*)、h₂(t) at the known angle η, the angle θ, the origin O, respectively_wCoordinate (w)_y,w_z) Said light source P_sCoordinate(s)_y,s_z) And said angleUnder the premise of (1), the projections of the steel balls at the corresponding positions t on the two bottom surfaces of the die body and the steel balls on the two bottom surfaces of the die bodyThe distance of the projection center of the nearest steel ball is tied down.

10. The method for calibrating a cone-beam CT imaging system as recited in claim 1,