CN113876346A - Iterative correction method for oblique image - Google Patents
Iterative correction method for oblique image Download PDFInfo
- Publication number
- CN113876346A CN113876346A CN202111352681.4A CN202111352681A CN113876346A CN 113876346 A CN113876346 A CN 113876346A CN 202111352681 A CN202111352681 A CN 202111352681A CN 113876346 A CN113876346 A CN 113876346A
- Authority
- CN
- China
- Prior art keywords
- panel detector
- flat panel
- inclination angle
- ellipse
- calculating
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 32
- 229910000831 Steel Inorganic materials 0.000 claims abstract description 66
- 239000010959 steel Substances 0.000 claims abstract description 66
- 230000011218 segmentation Effects 0.000 claims description 10
- 238000004364 calculation method Methods 0.000 claims description 8
- PHTXVQQRWJXYPP-UHFFFAOYSA-N ethyltrifluoromethylaminoindane Chemical compound C1=C(C(F)(F)F)C=C2CC(NCC)CC2=C1 PHTXVQQRWJXYPP-UHFFFAOYSA-N 0.000 claims description 4
- 238000001914 filtration Methods 0.000 claims description 3
- 238000010586 diagram Methods 0.000 description 6
- 239000000203 mixture Substances 0.000 description 2
- 239000007787 solid Substances 0.000 description 2
- 238000000605 extraction Methods 0.000 description 1
- 238000012163 sequencing technique Methods 0.000 description 1
Images
Classifications
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
- A61B6/02—Arrangements for diagnosis sequentially in different planes; Stereoscopic radiation diagnosis
- A61B6/03—Computed tomography [CT]
- A61B6/032—Transmission computed tomography [CT]
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
- A61B6/40—Arrangements for generating radiation specially adapted for radiation diagnosis
- A61B6/4064—Arrangements for generating radiation specially adapted for radiation diagnosis specially adapted for producing a particular type of beam
- A61B6/4085—Cone-beams
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
- A61B6/52—Devices using data or image processing specially adapted for radiation diagnosis
- A61B6/5211—Devices using data or image processing specially adapted for radiation diagnosis involving processing of medical diagnostic data
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
- A61B6/58—Testing, adjusting or calibrating thereof
- A61B6/582—Calibration
- A61B6/583—Calibration using calibration phantoms
Landscapes
- Health & Medical Sciences (AREA)
- Life Sciences & Earth Sciences (AREA)
- Engineering & Computer Science (AREA)
- Medical Informatics (AREA)
- Biomedical Technology (AREA)
- Heart & Thoracic Surgery (AREA)
- High Energy & Nuclear Physics (AREA)
- Physics & Mathematics (AREA)
- Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
- Optics & Photonics (AREA)
- Pathology (AREA)
- Radiology & Medical Imaging (AREA)
- Veterinary Medicine (AREA)
- Biophysics (AREA)
- Molecular Biology (AREA)
- Surgery (AREA)
- Animal Behavior & Ethology (AREA)
- General Health & Medical Sciences (AREA)
- Public Health (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Pulmonology (AREA)
- Theoretical Computer Science (AREA)
- Length-Measuring Devices Using Wave Or Particle Radiation (AREA)
- Analysing Materials By The Use Of Radiation (AREA)
Abstract
The invention discloses an iterative correction method of a tilted image, which comprises the following steps: (1) scanning the calibration die body by CT to obtain a projected image of the calibration die body on a detector, and acquiring the central coordinates of each steel ball in the projected image; (2) calculating parameters of an ellipse projected to the detector by two layers of steel balls; (3) calculating the current detector inclination angle according to the design parameters, the ellipse parameters and the distance from the ray source to the detector; (4) correcting the detector, and calculating the corrected elliptical parameters and the actual distance from the ray source to the detector; (5) calculating a corrected inclination angle of the detector; (6) judging whether the corrected inclination angle is smaller than a preset threshold value, if so, ending the circulation; and if not, taking the sum of the corrected inclination angle and the current inclination angle as the latest inclination angle, and returning to the step (4). According to the perspective model theory, the invention calculates the angle error by utilizing the change of the projection size caused by the rotation of the detector in different directions, thereby improving the robustness of angle solution while ensuring the precision.
Description
Technical Field
The invention relates to the technical field of image processing, in particular to an iterative correction method for a tilted image.
Background
The cone beam CT technology is a technology that an X-ray light source and a flat panel detector are adopted to carry out rotary scanning on an object to obtain projections of the object under different angles, and then a cone beam back projection reconstruction algorithm is used to obtain a three-dimensional image of the object. The implementation of the cone-beam back-projection reconstruction algorithm requires that it be provided with the geometrical parameters of the system consisting of the X-ray source, the detector and the object. The geometric parameters have important influence on the resolution, the artifacts and the like of the reconstructed image, and the geometric parameters need to be accurately calibrated for obtaining a high-quality reconstructed CT image.
The existing geometric parameter calibration method mainly depends on two-dimensional image information, and adopts a step-by-step solving geometric calibration method, and the calibration methods usually do not consider the problem of mutual influence among all geometric parameters, such as the influence of a detector under the condition of inclination in different directions and the mutual influence caused by deformation of CT equipment among different rotation angles; in addition, the prior calibration method mostly adopts a 'convergent point' to calculate the inclination angle, but the solving precision of the 'convergent point' fluctuates greatly due to the influence of the inclination of angles in different directions, so that the overall geometric parameter calibration precision is influenced.
Disclosure of Invention
The purpose of the invention is as follows: in order to solve the problems, the invention provides an iterative correction method of an oblique image, which improves the robustness of angle solution while ensuring the precision.
The technical scheme is as follows:
a method of iterative correction of a tilted image, comprising the steps of:
(1) scanning the calibration die body through CT to obtain a projection image of the calibration die body on a flat panel detector, and acquiring the central coordinates of steel balls in the image, wherein the steel balls are arranged in two layers and are arranged in the calibration die body in parallel, and each layer of steel balls form a standard circle;
(2) respectively calculating to obtain an equation of an ellipse formed on the projection image, and solving corresponding ellipse parameters according to the equation;
(3) calculating to obtain a current inclination angle of the current flat panel detector relative to an ideal position of the current flat panel detector according to the design parameters of the calibration die body, the elliptic parameters obtained in the step (2) and the distance from the ray source to the flat panel detector, wherein the initial value of the distance from the ray source to the flat panel detector is a design value;
(4) correcting the inclination angle of the image by using the current inclination angle of the flat panel detector, and calculating the corrected ellipse parameters and the actual distance from the ray source to the flat panel detector;
(5) calculating to obtain a corrected inclination angle of the flat panel detector relative to an ideal position after correction according to the design parameters of the calibration die body and the elliptic parameters obtained in the step (4) and the actual distance from the ray source to the flat panel detector;
(6) judging whether the inclination angle corrected in the step (5) is smaller than a preset threshold value, if so, ending the circulation; and (4) if not, taking the sum of the corrected inclination angle obtained in the step (5) and the current inclination angle in the step (4) as the latest current inclination angle, and returning to the step (4).
In the step (1), the number of steel balls on each layer in the calibration die body is set to be 4-16, and the diameter of each steel ball is set to be 1-5 mm.
In the calibration die body, the diameter of each steel ball is the same.
The number of steel balls in each layer is 12, and the diameter of each steel ball is 2 mm.
The calibration die body is a cylinder, and two layers of steel balls are symmetrically arranged on the upper side and the lower side of the center of the calibration die body.
In the step (1), the obtaining of the center coordinates of each steel ball specifically comprises:
the image is subjected to Gaussian filtering and threshold segmentation, a segmentation threshold is determined according to the number of pixel points projected by the steel balls, the position coordinates of the steel balls are preliminarily extracted through a self-adaptive threshold segmentation algorithm, the edges of the steel balls are extracted and segmented according to the position coordinates, the accurate outlines of the steel balls are obtained, and the central coordinates of the steel balls are obtained by performing ellipse fitting through a least square method.
The calculating of the inclination angle of the flat panel detector relative to the ideal position in the step (3) specifically comprises the following steps:
the inclination angle of the flat panel detector around the axis parallel to the vertical side of the flat panel detector is phi relative to the ideal position of the flat panel detector, and the inclination angle of the flat panel detector around the axis parallel to the horizontal side of the flat panel detector is theta;
respectively calculating according to the formulas (1) and (2) to obtain the corresponding inclination angle phiiAnd thetai；
Wherein, V0φ、V1φThe angle of inclination of the flat-panel detector relative to its ideal position about an axis parallel to its vertical edge is phiiUnder the condition of (1), the circle formed by the steel balls on the upper layer projects on the actual flat panel detector to form the vertex of the long axis of the ellipse; v2φ、V3φThe angle of inclination of the flat-panel detector relative to its ideal position about an axis parallel to its vertical edge is phiiUnder the condition of (1), the circle formed by the lower steel ball projects on an actual detector to form the vertex of the major axis of the ellipse; v0θ、V1θFor the inclination of the flat-panel detector with respect to its ideal position about an axis parallel to its horizontal sides at an angle thetaiUnder the condition of (3), the circle formed by the steel balls on the upper layer projects on an actual detector to form the vertex of the long axis of the ellipse; v2θ、V3θFor the inclination of the flat-panel detector with respect to its ideal position about an axis parallel to its horizontal sides at an angle thetaiUnder the condition of (1), the circle formed by the lower steel ball projects on an actual detector to form the vertex of the major axis of the ellipse; d2The initial value is the design value for the distance from the ray source to the flat panel detector.
The step (3) of calculating the inclination angle of the flat panel detector relative to the ideal position further comprises: at a corresponding inclination angle phiiAnd thetaiThe projection image is corrected in an inclined mode, the ellipse parameters in the corrected image are obtained, and the inclination angle eta of the flat panel detector around the ray source and taking the projection of the flat panel detector and the ray source connecting line as the axis is calculated according to the ellipse parametersi。
The angle of inclination etaiThe calculation of (a) is specifically:
V0φV1φand V2φV3φThe intersection point of the straight lines is Pφ，PφThe vertical point on the flat panel detector is PaThe point P is calculated by the formula (3)aDistance d to origin OOa：
Wherein d isO1Is the distance from the center of the upper ellipse to the origin O, dO2The distance from the center of the lower layer ellipse to the origin O is shown, a1 and b1 are coefficients of the upper layer ellipse equation, and a2 and b2 are coefficients of the lower layer ellipse equation;
calculating according to the formula (4) to obtain the inclination angle eta of the flat panel detector around the ray source, wherein the inclination angle eta is the axis of the projection of the flat panel detector on the flat panel detector and the connecting line of the ray sourcei：
η＝[d1aA(x，L1)+d2aA(x，L2)]/d12 (4)
Wherein L1 and L2 each represent V0φV1φAnd V2φV3φOn the straight line, A (x, L1) and A (x, L2) are respectively
And
the included angle between the straight line and the axis parallel to the horizontal edge of the flat panel detector is formed; d1a、d2aRespectively represents the distance from the center of the ellipse of the upper and lower layers to the point Pa, d12Indicating the distance between the centers of the upper and lower ellipses.
The calculation of the corrected ellipse corresponding parameters and the distance from the ray source to the actual flat panel detector in the step (5) is specifically as follows:
calculating a coordinate value of the ray source S in an actual coordinate system of the current flat panel detector according to the formula (5):
calculating the actual distance d2 from the ray source S to the actual flat panel detector according to the following formula:
the step (6) is specifically as follows:
calculating the actual distance d2 from the ray source S to the actual flat panel detector according to the current inclination angle in the step (4), and calculating the corrected inclination angle phi according to the ellipse parameter combination formulas (1), (2), (3) and (4) obtained in the step (4)i′、θi′、ηi′；
If phiiIf the absolute value of' is less than the first threshold, let phii' -0; if phiiIf the absolute value of' is not less than the first threshold value, let phii+1＝φi+φi' as the latest current tilt angle, and performing a tilt angle phi on the imagei' correction;
if thetai' if the absolute value is less than the first threshold, let θi' -0; if thetai' if the absolute value is not less than the first threshold, let θi+1＝θi+θi' as the latest current tilt angle, and subjecting the image to the tilt angle thetai' correction;
if etaiIf the absolute value of ` is less than the second threshold value, let ηi' -0; if etai' if the absolute value is not less than the second threshold, let ηi+1＝ηi+ηi' as the latest current tilt angle, and subjecting the image to a rotation angle ηi' of the above.
The first threshold is set to δ 0.00087rad, and the second threshold is set to δ 1 0.000087 rad.
And (5) when the corrected inclination angle is smaller than a preset threshold value, calculating a calibration parameter according to the current inclination angle, and outputting the calibration parameter.
Has the advantages that: according to the perspective model theory, the invention calculates the angle error by using the change of the projection size of the object caused by the rotation of the detector in different directions, thereby improving the robustness of angle solution while ensuring the precision. Meanwhile, in order to eliminate the problem of mutual influence among different parameters, the invention adopts an iterative calibration and successive correction mode to gradually correct the geometric parameter error in the calibration process until the requirement of calibration precision is met.
Drawings
FIG. 1 is a schematic projection diagram of an X-ray source, a calibration phantom and a flat panel detector;
FIG. 2 is a flow chart of the present invention;
FIG. 3 is a schematic view of a fixed mold body projection at a flat panel detector with a slope angle phi;
FIG. 4 is a schematic projection diagram of a fixed mold body when an inclination angle θ of the flat panel detector exists;
FIG. 5 is a schematic view of the flat panel detector fixed mold body projection with the inclination angle phi and theta simultaneously;
FIG. 6 is a schematic diagram of the module projection when the flat panel detector has an inclination angle eta.
Detailed Description
The invention is further elucidated with reference to the drawings and the embodiments.
In the embodiment disclosed by the invention, the design parameters of the calibration die body are as follows: the calibration die body is cylindrical, two layers of steel balls are arranged in the calibration die body in parallel, the number of each layer of steel balls is 12, the diameter of each steel ball is 2mm, 12 steel balls in each layer are arranged on a standard circular track, the standard circular track is symmetrically arranged on the upper side and the lower side of the center of the calibration die body, so that the two layers of steel balls are symmetrically arranged on the upper side and the lower side of the center of the calibration die body, and the radius of the standard circular track is R and the unit of mm; the vertical distance between the upper and lower standard circular tracks is D in mm. The invention adopts an X-ray source of cone-beam CT to scan a calibration die body and project the calibration die body to a flat panel detector, the X-ray source penetrates through the calibration die body, the central line of rays emitted by the X-ray source is vertical to the flat panel detector, but in the actual operation process, the flat panel detector has angle deviation compared with the ideal position, and the invention obtains the rotation angle of the flat panel detector compared with the ideal position through iterative calculation so as to correct the geometric parameter error.
Specifically, before correction, the calibration phantom is placed in the center of the field of view, the calibration phantom is adjusted to the center of the field of view by taking images at the right and side positions, then the calibration phantom is rotationally scanned by cone beam CT, the rotation axis of the cone beam CT is basically coincident with the center line of the calibration phantom to obtain a sequence image, as shown in fig. 1, an X-ray beam emitted by an X-ray source S passes through the calibration phantom, and the X-ray beam is projected on a flat panel detector to generate a sequence imageAnd projecting the image. The projection of a standard circle formed by the upper and lower layers of steel balls on the plane detector is an ellipse, the top points of the upper and lower ellipses which are respectively farthest away, namely the top points of the major axes of the ellipses are respectively V0、V1、V2、V3。
Fig. 2 is a flowchart of an iterative correction method for a tilted image according to the present invention, and as shown in fig. 2, the iterative correction method for a tilted image according to the present invention includes the following steps:
(1) carrying out rotary scanning on the calibration die body through the cone beam CT to obtain a projection image of the calibration die body on the flat panel detector;
(2) acquiring the center coordinates of each steel ball: performing Gaussian filtering and threshold segmentation on a single projected image, determining a segmentation threshold according to the number of pixel points projected by the steel balls, preliminarily extracting the position coordinates of the steel balls through a self-adaptive threshold segmentation algorithm, performing edge extraction and segmentation on the steel balls according to the position coordinates, acquiring the accurate outline of each steel ball, and performing ellipse fitting by adopting a least square method to acquire the central coordinates of each steel ball;
(3) sequencing the steel balls on the upper layer and the lower layer obtained in the step (2) in a counterclockwise manner, and performing centrosymmetric pairing on the steel balls on the upper layer and the lower layer, wherein the intersection point between connecting lines of each pair of centrosymmetric steel balls is the center of the calibration phantom, and the coordinate of the central point of the calibration phantom projected on the flat panel detector is obtained by solving according to the center point;
(4) respectively performing least square ellipse fitting on the upper and lower steel ball sets to respectively obtain an ellipse equation of the upper and lower steel balls, and solving to obtain a long axis vertex and an ellipse center coordinate of an ellipse formed by projection of the upper and lower steel balls on the flat panel detector;
(5) calculating the inclination angle of the flat panel detector;
FIG. 3 is a schematic view of a projection of a calibration phantom onto a flat panel detector when the detector has an inclination angle φ, as shown in FIG. 3, a dashed rectangle represents an ideal position of the flat panel detector, and a solid rectangle represents an actual position of the flat panel detector; when the flat panel detector is in an ideal position, the projection of a ray of the X-ray source S passing through the center of the calibration die body on the flat panel detector is O; establishing an ideal coordinate system by taking the O point as an origin, taking an axis which passes through the O point and is parallel to the vertical edge of the flat panel detector as a y axis, taking an axis which passes through the origin O and is parallel to the horizontal edge of the flat panel detector as an x axis, and taking an axis which passes through the O point and is simultaneously perpendicular to the x axis and the y axis as a z axis;
when the inclination angle of the flat panel detector relative to the ideal position of the flat panel detector around the y-axis is phi, an included angle exists between the actual coordinate system where the flat panel detector is located and the ideal coordinate system, and at this time, the actual coordinate system is formed by xφAxis, yφAxis and zφShaft composition, wherein xφThe axis is the x axis and is obtained by rotating an angle phi around the corresponding direction of the y axisφThe axis being coincident with the y-axis, zφRotating the axis which is the z axis around the corresponding direction of the y axis by a rotation angle phi;
fig. 4 is a schematic projection diagram of a calibration mold body on a flat panel detector when an inclination angle θ exists in the flat panel detector, as shown in fig. 4, a dashed rectangle represents an ideal position of the flat panel detector under an ideal condition, and a solid rectangle represents an actual position of the flat panel detector; when the inclination angle theta of the flat panel detector around the x axis is relative to the ideal position of the flat panel detector, an included angle exists between the actual coordinate system of the flat panel detector and the ideal coordinate system of the flat panel detector, and the actual coordinate system is formed by xθAxis, yθAxis and zθShaft composition, wherein xθThe axis being coincident with the x-axis, yθRotation angle theta around x axis with axis y as well as zθThe axis is the rotation angle theta of the z axis around the x axis;
FIG. 5 is a schematic diagram of a fixed mold body projection when the detector has inclination angle phi and theta at the same time, and the current inclination angle values phi and theta of the flat panel detector are obtained by calculation according to the formulas (1) and (2)iAnd thetai；
Wherein, V0φ、V1φFor flat panel probingThe angle of inclination of the device about the y-axis with respect to the ideal position of the flat panel detector is phiiUnder the condition of (1), the circle formed by the steel balls on the upper layer projects on the actual flat panel detector to form the vertex of the long axis of the ellipse; v2φ、V3φThe inclination angle around the y-axis for the ideal position of the flat panel detector relative to the flat panel detector is phiiUnder the condition of (1), the circle formed by the lower steel ball projects on the actual flat panel detector to form the vertex of the major axis of the ellipse; v0θ、V1θThe inclination angle around the x axis of the flat panel detector relative to the ideal position of the flat panel detector is thetaiUnder the condition of (1), the circle formed by the steel balls on the upper layer projects on the actual flat panel detector to form the vertex of the long axis of the ellipse; v2θ、V3θThe inclination angle around the x axis of the flat panel detector relative to the ideal position of the flat panel detector is thetaiUnder the condition of (1), the circle formed by the lower steel ball projects on the actual flat panel detector to form the vertex of the major axis of the ellipse; d1Is the distance from the X-ray source S to the center of the calibration phantom, d2The distance from the ray source S to the flat panel detector is calculated for the first time, d1、d2For design values, during later iterations, d1、d2The value of (d) is the actual distance calculated in equation (7);
the current flat panel detector existing inclination angle phi obtained by the calculationiAnd thetaiThe projection image is subjected to tilt correction, ellipse parameters in the corrected image are obtained, and a current value eta of a tilt angle eta of the flat panel detector relative to the ideal position of the flat panel detector around the z axis is calculated according to the ellipse parametersiFig. 6 is a schematic diagram of a projection of a fixed mold body when the flat panel detector has an η -scale inclination angle, which is specifically calculated as follows:
as shown in FIG. 3, V0φV1φAnd V2φV3φThe intersection point of the straight lines is Pφ，PφAt yφThe vertical point on the shaft is PaIt can be calculated from equation (3):
wherein，dOaIs a point PaDistance to origin O, dO1Is the distance from the center of the upper ellipse to the origin O, dO2Is the distance from the center of the lower layer ellipse to the origin O, a1、b1Is the coefficient of the upper elliptic equation, a2、b2The coefficients of the lower elliptic equation;
calculating a current value eta of a tilt angle eta of the ideal position of the flat panel detector relative to the flat panel detector around the z-axis according to the formula (4)i；
η＝[d1aA(x，L1)+d2aA(x，L2)]/d12 (4)
Wherein L is1、L2Respectively represent V0φV1φAnd V2φV3φThe straight line, A (x, L)1) And A (x, L)2) Are respectively as
And
the included angle between the straight line and the x axis; d1a、d2aRespectively represents the center to point P of the ellipse of the upper and lower layersaDistance of d12Represents the distance between the centers of the upper and lower ellipses;
(6) calculating the actual distance d from the X-ray source S to the center of the calibration phantom after the flat panel detector is corrected1And the actual distance d from the X-ray source S to the flat panel detector2；
When there is a tilt angle value phi according to the ideal position of the flat panel detector relative to the flat panel detectori、θiAnd ηiCalculating the coordinate values of the rotation center Q (namely the center of the calibration phantom), the world coordinate system of the X-ray source S and the actual coordinate system of the flat panel detector, and calculating the actual distance d from the X-ray source S to the center of the calibration phantom according to the coordinate values1And the actual distance d from the X-ray source S to the flat panel detector2(ii) a The method specifically comprises the following steps:
(61) calculating the coordinate value of the X-ray source S in the actual coordinate system of the current flat panel detector according to the formula (5):
and meanwhile, calculating a coordinate value of the rotation center Q in an actual coordinate system of the current flat panel detector according to the formula (6):
(62) calculating according to the formula (7) to obtain the actual distance d from the X-ray source S to the center of the calibration phantom1And the actual distance d from the X-ray source S to the flat panel detector2：
(7) Calculating a current value phii、θi、ηiCorrected inclination angle phi between corrected flat panel detector and ideal position thereofi′、θi、ηi', and comparing with a set threshold;
calculating the actual distance d from the new X-ray source S to the center of the calibration phantom according to the formula (7)1And the actual distance d from the X-ray source S to the flat panel detector2Respectively calculating the corrected inclination angle phi of the corrected flat panel detector by combining the formulas (1) and (2)i′、θi′；
Setting the first threshold δ to 0.00087rad, if φiIf the absolute value of' is less than the first threshold, let phii' -0; if phiiIf the absolute value of' is not less than the first threshold value, let phii+1＝φi+φi' and performing inclination angle phi on the imagei' correction is performed;
if thetai' if the absolute value is less than the first threshold, let θi' -0; if thetai' if the absolute value is not less than the first threshold, let θi+1＝θi+θi', and to the image tilt angle thetai' correction is performed;
calculating the inclination angle eta of the corrected flat panel detector by the formula (4)i' setting a second threshold value delta10.000087rad, if ηiIf the absolute value of ` is less than the second threshold value, let ηi' -0; if etai' if the absolute value is not less than the second threshold, let ηi+1＝ηi+ηi', and to the image tilt angle etai' correction is performed;
(8) judgment of phii′、θi′、ηiIf yes, ending the circulation, taking the current inclination angle as an included angle between the actual flat panel detector and the ideal detector position to calculate a calibration parameter and output the calibration parameter, wherein the calibration parameter is a rotation center obtained by calculation according to the current inclination angle, a coordinate value of the ray source in a flat panel detector coordinate system, and an actual distance d from the ray source to the center of a calibration die body1And the actual distance d from the ray source to the flat panel detector2(ii) a Otherwise, with phii+1、θi+1、ηi+1As the latest current value phii、θi、ηiAnd (6) returning.
In other embodiments, the number of the steel balls on each layer can be set to be 4-16, the diameter of the steel balls can be set to be 1-5 mm, and the centers of the steel balls are on the same circumference; preferably, all the steel balls are equal in diameter.
In the present invention, the specific values of the first threshold and the second threshold may be set according to actual conditions.
According to the perspective model theory, the invention calculates the angle error by utilizing the change of the projection size of the object caused by the rotation of the flat panel detector in different directions, thereby improving the robustness of angle solution while ensuring the precision. Meanwhile, in order to eliminate the problem of mutual images among different parameters, the invention adopts the modes of iterative calibration and successive correction to gradually correct the geometric parameter errors until the requirement of calibration precision is met.
Although the preferred embodiments of the present invention have been described in detail, the present invention is not limited to the details of the foregoing embodiments, and various equivalent changes (such as number, shape, position, etc.) may be made to the technical solution of the present invention within the technical spirit of the present invention, and these equivalent changes are all within the protection scope of the present invention.
Claims (13)
1. A method of iterative correction of a tilted image, characterized by: the method comprises the following steps:
(1) scanning the calibration die body through CT to obtain a projected image of the calibration die body on a flat panel detector, and acquiring the central coordinates of steel balls in the projected image, wherein the steel balls are arranged in two layers and are arranged in the calibration die body in parallel, and each layer of steel balls form a standard circle;
(2) respectively calculating to obtain an equation of an ellipse formed on the projection image, and solving corresponding ellipse parameters according to the equation;
(3) calculating to obtain a current inclination angle of the current flat panel detector relative to an ideal position of the current flat panel detector according to the design parameters of the calibration die body, the elliptic parameters obtained in the step (2) and the distance from the ray source to the flat panel detector, wherein the initial value of the distance from the ray source to the flat panel detector is a design value;
(4) correcting the inclination angle of the image by using the current inclination angle of the flat panel detector, and calculating the corrected ellipse parameters and the actual distance from the ray source to the flat panel detector;
(5) calculating to obtain a corrected inclination angle of the flat panel detector relative to an ideal position after correction according to the design parameters of the calibration die body and the elliptic parameters obtained in the step (4) and the actual distance from the ray source to the flat panel detector;
(6) judging whether the inclination angle corrected in the step (5) is smaller than a preset threshold value, if so, ending the circulation; and (4) if not, taking the sum of the corrected inclination angle obtained in the step (5) and the current inclination angle in the step (4) as the latest current inclination angle, and returning to the step (4).
2. The iterative correction method for a tilted image according to claim 1, characterized in that: in the step (1), the number of steel balls on each layer in the calibration die body is set to be 4-16, and the diameter of each steel ball is set to be 1-5 mm.
3. The iterative correction method for a tilted image according to claim 2, characterized in that: in the calibration die body, the diameter of each steel ball is the same.
4. The iterative correction method for a tilted image according to claim 3, characterized in that: the number of steel balls in each layer is 12, and the diameter of each steel ball is 2 mm.
5. The iterative correction method for a tilted image according to claim 1, characterized in that: the calibration die body is a cylinder, and two layers of steel balls are symmetrically arranged on the upper side and the lower side of the center of the calibration die body.
6. The iterative correction method for a tilted image according to claim 1, characterized in that: in the step (1), the obtaining of the center coordinates of each steel ball specifically comprises:
the image is subjected to Gaussian filtering and threshold segmentation, a segmentation threshold is determined according to the number of pixel points projected by the steel balls, the position coordinates of the steel balls are preliminarily extracted through a self-adaptive threshold segmentation algorithm, the edges of the steel balls are extracted and segmented according to the position coordinates, the accurate outlines of the steel balls are obtained, and the central coordinates of the steel balls are obtained by performing ellipse fitting through a least square method.
7. The iterative correction method for a tilted image according to claim 1, characterized in that: the calculating of the inclination angle of the flat panel detector relative to the ideal position in the step (3) specifically comprises the following steps:
the inclination angle of the flat panel detector around the axis parallel to the vertical side of the flat panel detector is phi relative to the ideal position of the flat panel detector, and the inclination angle of the flat panel detector around the axis parallel to the horizontal side of the flat panel detector is theta;
respectively calculating according to the formulas (1) and (2) to obtain the corresponding inclination angle phiiAnd thetai；
Wherein, V0φ、V1φThe angle of inclination of the flat-panel detector relative to its ideal position about an axis parallel to its vertical edge is phiiUnder the condition of (1), the circle formed by the steel balls on the upper layer projects on the actual flat panel detector to form the vertex of the long axis of the ellipse; v2φ、V3φThe angle of inclination of the flat-panel detector relative to its ideal position about an axis parallel to its vertical edge is phiiUnder the condition of (1), the circle formed by the lower steel ball projects on an actual detector to form the vertex of the major axis of the ellipse; v0θ、V1θFor the inclination of the flat-panel detector with respect to its ideal position about an axis parallel to its horizontal sides at an angle thetaiUnder the condition of (3), the circle formed by the steel balls on the upper layer projects on an actual detector to form the vertex of the long axis of the ellipse; v2θ、V3θFor the inclination of the flat-panel detector with respect to its ideal position about an axis parallel to its horizontal sides at an angle thetaiUnder the condition of (1), the circle formed by the lower steel ball projects on an actual detector to form the vertex of the major axis of the ellipse; d2The initial value is the design value for the distance from the ray source to the flat panel detector.
8. The iterative correction method for a tilted image according to claim 7, characterized in that: the step (3) of calculating the inclination angle of the flat panel detector relative to the ideal position further comprises: at a corresponding inclination angle phiiAnd thetaiThe projection image is corrected in an inclined mode, the ellipse parameters in the corrected image are obtained, and the inclination angle eta of the flat panel detector around the ray source and taking the projection of the flat panel detector and the ray source connecting line as the axis is calculated according to the ellipse parametersi。
9. The iterative correction method for a tilted image according to claim 8, characterized in that: the angle of inclination etaiThe calculation of (a) is specifically:
V0φV1φand V2φV3φThe intersection point of the straight lines is Pφ，PφThe vertical point on the flat panel detector is PaThe point P is calculated by the formula (3)aDistance d to origin OOa：
Wherein d isO1Is the distance from the center of the upper ellipse to the origin O, dO2Is the distance from the center of the lower layer ellipse to the origin O, a1、b1Is the coefficient of the upper elliptic equation, a2、b2The coefficients of the lower elliptic equation;
calculating according to the formula (4) to obtain the inclination angle eta of the flat panel detector around the ray source, wherein the inclination angle eta is the axis of the projection of the flat panel detector on the flat panel detector and the connecting line of the ray sourcei：
η＝[d1aA(x，L1)+d2aA(x，L2)]/d12 (4)
Wherein L is1、L2Respectively represent V0φV1φAnd V2φV3φOn the straight line, A (x, L)1) And A (x, L)2) Are respectively as
And
the included angle between the straight line and the axis parallel to the horizontal edge of the flat panel detector is formed; d1a、d2aRespectively represents the center to point P of the ellipse of the upper and lower layersaDistance of d12Indicating the distance between the centers of the upper and lower ellipses.
10. The iterative correction method for a tilted image according to claim 9, characterized in that: the calculation of the corrected ellipse corresponding parameters and the distance from the ray source to the actual flat panel detector in the step (5) is specifically as follows:
calculating a coordinate value of the ray source S in an actual coordinate system of the current flat panel detector according to the formula (5):
calculating the actual distance d from the ray source S to the actual flat panel detector according to the following formula2：
11. The iterative correction method for a tilted image according to claim 10, characterized in that: the step (6) is specifically as follows:
calculating the actual distance d from the ray source S to the actual flat panel detector according to the current inclination angle in the step (4)2And calculating to obtain a corrected inclination angle phi according to the ellipse parameters obtained in the step (4) and the combination formulas (1), (2), (3) and (4)i′、θi′、ηi′；
If phiiIf the absolute value of' is less than the first threshold, let phii' -0; if phiiIf the absolute value of' is not less than the first threshold value, let phii+1＝φi+φi' as the latest current tilt angle, and performing a tilt angle phi on the imagei' correction;
if thetai' if the absolute value is less than the first threshold, let θi' -0; if thetai' if the absolute value is not less than the first threshold, let θi+1＝θi+θi' as the latest current tilt angle, and subjecting the image to the tilt angle thetai' correction;
if etaiIf the absolute value of ` is less than the second threshold value, let ηi' -0; if etai' if the absolute value is not less than the second threshold, let ηi+1＝ηi+ηi' as the latest current tilt angle, and rotates the angle of the imageηi' of the above.
12. The iterative correction method for a tilted image according to claim 10, characterized in that: the first threshold is set to be δ 0.00087rad, and the second threshold is set to be δ1＝0.000087rad。
13. The iterative correction method for a tilted image according to any one of claims 1 to 12, characterized in that: and (5) when the corrected inclination angle is smaller than a preset threshold value, calculating a calibration parameter according to the current inclination angle, and outputting the calibration parameter.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111352681.4A CN113876346B (en) | 2021-11-16 | 2021-11-16 | Iterative correction method for inclined image |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111352681.4A CN113876346B (en) | 2021-11-16 | 2021-11-16 | Iterative correction method for inclined image |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN113876346A true CN113876346A (en) | 2022-01-04 |
| CN113876346B CN113876346B (en) | 2024-06-07 |
Family
ID=79017527
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202111352681.4A Active CN113876346B (en) | 2021-11-16 | 2021-11-16 | Iterative correction method for inclined image |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN113876346B (en) |
Cited By (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115105757A (en) * | 2022-05-24 | 2022-09-27 | 沈阳东软智睿放疗技术有限公司 | Medical die body, and equipment precision detection method and system based on medical die body |
| CN115963130A (en) * | 2022-12-29 | 2023-04-14 | 中国科学院福建物质结构研究所 | A Correction Method for Non-Gaussian Deviation of High Angle Diffraction Data in X-ray Single Crystal Diffraction Experiment |
| CN118172397A (en) * | 2024-05-09 | 2024-06-11 | 杭州脉流科技有限公司 | Three-dimensional inclination correction method, computer equipment and computer program product for intracranial CT (computed tomography) image based on image registration |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US5293312A (en) * | 1991-06-20 | 1994-03-08 | Waggener Robert G | Method and apparatus for computing tomographic scans |
| US20070172025A1 (en) * | 2006-01-26 | 2007-07-26 | Kabushiki Kaisha Toshiba | X-ray ct apparatus, method of aligning phantom, and phantom retaining tool |
| US20150216498A1 (en) * | 2012-08-20 | 2015-08-06 | Orangedental Gmbh & Co. Kg | Geometric Characterization and Calibration of a Cone-Beam Computer Tomography Apparatus |
| US20170020481A1 (en) * | 2013-11-28 | 2017-01-26 | Nikon Metrology Nv | Calibration apparatus and method for computed tomography |
| CN108122203A (en) * | 2016-11-29 | 2018-06-05 | 上海东软医疗科技有限公司 | A kind of bearing calibration of geometric parameter, device, equipment and system |
-
2021
- 2021-11-16 CN CN202111352681.4A patent/CN113876346B/en active Active
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US5293312A (en) * | 1991-06-20 | 1994-03-08 | Waggener Robert G | Method and apparatus for computing tomographic scans |
| US20070172025A1 (en) * | 2006-01-26 | 2007-07-26 | Kabushiki Kaisha Toshiba | X-ray ct apparatus, method of aligning phantom, and phantom retaining tool |
| US20150216498A1 (en) * | 2012-08-20 | 2015-08-06 | Orangedental Gmbh & Co. Kg | Geometric Characterization and Calibration of a Cone-Beam Computer Tomography Apparatus |
| US20170020481A1 (en) * | 2013-11-28 | 2017-01-26 | Nikon Metrology Nv | Calibration apparatus and method for computed tomography |
| CN108122203A (en) * | 2016-11-29 | 2018-06-05 | 上海东软医疗科技有限公司 | A kind of bearing calibration of geometric parameter, device, equipment and system |
Non-Patent Citations (1)
| Title |
|---|
| 张震: "锥束口腔CT 机设计及其关键技术研究", 工程科技Ⅱ辑, no. 06, pages 1 - 71 * |
Cited By (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115105757A (en) * | 2022-05-24 | 2022-09-27 | 沈阳东软智睿放疗技术有限公司 | Medical die body, and equipment precision detection method and system based on medical die body |
| CN115105757B (en) * | 2022-05-24 | 2025-07-18 | 沈阳东软智睿放疗技术有限公司 | Medical die body, and equipment precision detection method and system based on medical die body |
| CN115963130A (en) * | 2022-12-29 | 2023-04-14 | 中国科学院福建物质结构研究所 | A Correction Method for Non-Gaussian Deviation of High Angle Diffraction Data in X-ray Single Crystal Diffraction Experiment |
| CN115963130B (en) * | 2022-12-29 | 2023-09-29 | 中国科学院福建物质结构研究所 | Correction method for non-Gaussian deviation of high-angle diffraction data of X-ray single crystal diffraction experiment |
| CN118172397A (en) * | 2024-05-09 | 2024-06-11 | 杭州脉流科技有限公司 | Three-dimensional inclination correction method, computer equipment and computer program product for intracranial CT (computed tomography) image based on image registration |
| CN118172397B (en) * | 2024-05-09 | 2024-07-26 | 杭州脉流科技有限公司 | Three-dimensional inclination correction method, computer equipment and computer program product for intracranial CT (computed tomography) image based on image registration |
Also Published As
| Publication number | Publication date |
|---|---|
| CN113876346B (en) | 2024-06-07 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN113876346A (en) | Iterative correction method for oblique image | |
| CN108122203B (en) | Geometric parameter correction method, device, equipment and system | |
| CN106667512B (en) | Geometric correction method of X-ray imaging device and breast tomography device | |
| WO2017181471A1 (en) | Calibration method and system for geometric calibration phantom | |
| CN109498051B (en) | Automatic position calibration method and system for CT sickbed frame | |
| CN102652674B (en) | A method and system for removing geometric artifacts in CT images | |
| WO2018218611A1 (en) | Geometric parameter determination method for cone beam computed tomography system | |
| CN102915535A (en) | Method and system for correcting circle center deviation of round mark points during camera projection transformation | |
| CN106706675B (en) | A Calibration Method Based on Computer Layered Scanning Imaging CL System | |
| CN117392211B (en) | BGA element rapid identification positioning method and system and storage medium | |
| CN108460738A (en) | Medical image sloped correcting method based on B-spline | |
| CN112603346A (en) | Detector deflection correction method based on marker imaging | |
| CN110084855A (en) | A kind of improvement CBCT geometrical parameter calibration algorithm | |
| CN115105757B (en) | Medical die body, and equipment precision detection method and system based on medical die body | |
| CN115836875A (en) | Calibration method and system | |
| CN112971984B (en) | Coordinate registration method based on integrated surgical robot | |
| CN110645928A (en) | Space coordinate positioning method of three-dimensional scanner | |
| CN118014911A (en) | Correction phantom, CBCT geometric correction method, system and device based on offset detector | |
| CN117137510A (en) | Correction method and device for geometrical parameter errors of cone beam CT system | |
| CN113643428A (en) | A full-parameter geometric calibration method suitable for multi-degree-of-freedom cone-beam CT | |
| CN118037605A (en) | Planar imaging perspective deformation correction method based on extrinsic parameter adjustment | |
| CN203763103U (en) | Geometric correction device of detector of cone beam CT (computed Tomography) system | |
| CN116681600A (en) | CBCT calibration method aiming at non-ideal circular track and poor in motion repeatability | |
| CN110501360B (en) | Standard device for correcting pose of micro CT system and implementation method | |
| CN117974826B (en) | Cone beam CT geometric parameter calibration method, device, equipment and storage medium |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| CB02 | Change of applicant information |
Address after: 210000 building 3, No. 34, Dazhou Road, Yuhuatai District, Nanjing, Jiangsu Province Applicant after: Tuodao Medical Technology Co.,Ltd. Address before: Room 102-86, building 6, 57 Andemen street, Yuhuatai District, Nanjing, Jiangsu 210000 Applicant before: Nanjing Tuodao Medical Technology Co.,Ltd. |
|
| CB02 | Change of applicant information | ||
| GR01 | Patent grant |