CN113876346B - Iterative correction method for inclined image - Google Patents
Iterative correction method for inclined image Download PDFInfo
- Publication number
- CN113876346B CN113876346B CN202111352681.4A CN202111352681A CN113876346B CN 113876346 B CN113876346 B CN 113876346B CN 202111352681 A CN202111352681 A CN 202111352681A CN 113876346 B CN113876346 B CN 113876346B
- Authority
- CN
- China
- Prior art keywords
- flat panel
- panel detector
- inclination angle
- ellipse
- detector
- 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.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 30
- 229910000831 Steel Inorganic materials 0.000 claims abstract description 71
- 239000010959 steel Substances 0.000 claims abstract description 71
- 230000011218 segmentation Effects 0.000 claims description 12
- 230000005855 radiation Effects 0.000 claims description 6
- 238000004364 calculation method Methods 0.000 claims description 4
- 238000000605 extraction Methods 0.000 claims description 3
- 238000001914 filtration Methods 0.000 claims description 3
- 238000002591 computed tomography Methods 0.000 abstract 2
- 238000010586 diagram Methods 0.000 description 3
- 239000007787 solid Substances 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000012163 sequencing technique Methods 0.000 description 1
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 an inclined image, which comprises the following steps: (1) The CT scanning calibration die body obtains a projection image of the CT scanning calibration die body on the detector, and the center coordinates of all steel balls in the projection image are obtained; (2) Calculating parameters of ellipses projected onto the detector by the 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 between the ray source and the detector; (5) calculating a detector corrected caster angle; (6) Judging whether the corrected back dip angle is smaller than a preset threshold value, if so, ending the cycle; otherwise, the sum of the corrected back tilt angle and the current tilt angle is taken as the latest tilt angle, and the step (4) is returned. According to the perspective model theory, the angle error is calculated by utilizing the change of the projection size caused by the rotation of the detector in different directions, so that the accuracy is ensured, and the robustness of angle solving is improved.
Description
Technical Field
The invention relates to the technical field of image processing, in particular to an iterative correction method of an inclined image.
Background
The cone beam CT technology is a technology for rotationally scanning an object by adopting an X-ray light source and a flat panel detector to obtain projections of the object under different angles, and then obtaining a three-dimensional image of the object by adopting a cone beam back projection reconstruction algorithm. Implementing a cone beam back projection reconstruction algorithm requires providing it with the geometrical parameters of the system consisting of the X-ray source, detector and object. These geometric parameters have important effects on resolution and artifacts of the reconstructed image, and for obtaining a high quality reconstructed CT image, precise calibration of these geometric parameters is required.
The conventional geometric parameter calibration method mainly depends on two-dimensional image information, and a step-by-step solution geometric calibration method is adopted, so that the mutual influence problem among various geometric parameters, such as the influence of the detector under the condition of inclination in different directions and the mutual influence caused by deformation of CT equipment among different rotation angles, is not considered in the calibration methods; in addition, the existing calibration method mostly adopts a converging point to calculate the inclination angle, but the solving precision of the converging point has larger fluctuation due to the influence of the inclination angles in different directions, thereby influencing the overall geometric parameter calibration precision.
Disclosure of Invention
The invention aims to: in order to solve the problems, the invention provides an iterative correction method of an inclined image, which improves the robustness of angle solving while ensuring the accuracy.
The technical scheme is as follows:
an iterative correction method of an oblique image, comprising the steps of:
(1) Scanning the calibration die body through CT to obtain a projection image of the calibration die body on the flat panel detector, and obtaining the center coordinates of steel balls in the image, wherein the steel balls are arranged in two layers and are arranged in parallel in the calibration die body, and each layer of steel balls form a standard circle;
(2) Respectively calculating equations of ellipses formed on the projection images, and solving corresponding ellipse parameters according to the equations;
(3) Calculating 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 ellipse parameters obtained in the step (2) and the distance between the radiation source and the flat panel detector, wherein an initial value of the distance between the radiation source and 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 corrected ellipse parameters and the actual distance from the ray source to the flat panel detector;
(5) Obtaining ellipse parameters and actual distance between the ray source and the flat panel detector according to the design parameters of the calibration die body and the step (4), and calculating to obtain corrected inclination angles of the flat panel detector relative to the ideal positions of the flat panel detector after correction;
(6) Judging whether the inclination angle corrected in the step (5) is smaller than a preset threshold value, if so, ending the cycle; otherwise, the sum of the corrected inclination angle obtained in the step (5) and the current inclination angle in the step (4) is used as the latest current inclination angle, and the step (4) is returned.
In the step (1), the number of steel balls in each layer in the calibration die body is set to be 4-16, and the diameters of the steel balls are 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 set to 12, and the diameter of each steel ball is set to 2mm.
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 center coordinates of the steel balls are obtained specifically as follows:
And carrying out Gaussian filtering and threshold segmentation on the image, determining a segmentation threshold value according to the number of pixel points projected by the steel balls, initially extracting the position coordinates of the steel balls through a self-adaptive threshold segmentation algorithm, carrying out edge extraction segmentation on the steel balls according to the position coordinates, obtaining the accurate contour of each steel ball, and carrying out ellipse fitting by adopting a least square method to obtain the center coordinates of each steel ball.
The inclination angle of the flat panel detector relative to the ideal position in the step (3) is specifically:
The flat panel detector has a tilt angle phi about an axis parallel to its vertical side and a tilt angle theta about an axis parallel to its horizontal side with respect to its ideal position;
Respectively calculating corresponding dip angles phi i and theta i according to formulas (1) and (2);
Wherein V 0φ、V1φ is the vertex of the long axis of the ellipse formed by the projection of the circle formed by the upper layer steel ball on the actual flat panel detector under the condition that the inclination angle of the flat panel detector relative to the ideal position of the flat panel detector around the axis parallel to the vertical edge of the flat panel detector is phi i; v 2φ、V3φ is the vertex of the long axis of the ellipse formed by the projection of the circle formed by the lower steel ball on the actual detector under the condition that the inclination angle of the flat panel detector relative to the ideal position of the flat panel detector around the axis parallel to the vertical edge of the flat panel detector is phi i; v 0θ、V1θ is the vertex of the long axis of the ellipse formed by the projection of the circle formed by the upper layer steel ball on the actual detector under the condition that the inclination angle of the flat panel detector relative to the ideal position of the flat panel detector around the axis parallel to the horizontal edge of the flat panel detector is theta i; v 2θ、V3θ is the vertex of the long axis of the ellipse formed by the projection of the circle formed by the lower steel ball on the actual detector under the condition that the inclination angle of the flat panel detector relative to the ideal position of the flat panel detector around the axis parallel to the horizontal edge of the flat panel detector is theta i; d 2 is the distance from the radiation source to the flat panel detector, and its initial value is the design value.
The calculating of the tilt angle of the flat panel detector relative to its ideal position in the step (3) further comprises: and (3) performing inclination correction on the projection image by using corresponding inclination angles phi i and theta i to obtain ellipse parameters in the corrected image, and calculating the inclination angle eta i of the flat panel detector diffraction line source taking the connection line of the flat panel detector projection and the ray source as an axis according to the ellipse parameters.
The calculation of the inclination angle eta i specifically comprises the following steps:
The intersection point of the straight lines where V 0φV1φ and V 2φV3φ are located is P φ,Pφ, the vertical point of the flat panel detector is P a, and the distance d Oa from the point P a to the origin O is calculated by the formula (3):
Wherein d O1 is the distance from the center of the upper ellipse to the origin O, d O2 is the distance from the center of the lower ellipse to the origin O, a1 and b1 are coefficients of the upper ellipse equation, and a2 and b2 are coefficients of the lower ellipse equation;
Calculating according to the formula (4), obtaining the inclination angle eta i of the flat panel detector diffraction line source taking the projection of the flat panel detector and the line connecting with the ray source as the axis:
η＝[d1aA(x,L1)+d2aA(x,L2)]/d12 (4)
wherein L1 and L2 respectively represent straight lines where V 0φV1φ and V 2φV3φ are located, and A (x, L1) and A (x, L2) are respectively And/>The included angle between the straight line and the axis parallel to the horizontal side of the flat panel detector; d 1a、d2a represents the distance from the center of the upper and lower ellipse to the point Pa, and d 12 represents the distance between the centers of the upper and lower ellipse.
The step (5) is to calculate the corrected ellipse corresponding parameter and the distance between the ray source and the actual flat panel detector specifically as follows:
calculating the coordinate value of the ray source S in the actual coordinate system where the current flat panel detector is located according to the formula (5):
the actual distance d2 of the source S to the actual flat panel detector is calculated according to the following equation:
The step (6) specifically comprises the following steps:
Calculating the actual distance d2 from the radiation 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 i′、θi′、ηi' according to the ellipse parameters obtained in the step (4) and the combination formulas (1), (2) (3) and (4);
If the absolute value of phi i 'is less than the first threshold, let phi i' =0; if the absolute value of phi i ' is not less than the first threshold value, phi i+1＝φi+φi ' is taken as the latest current inclination angle, and the inclination angle phi i ' of the image is corrected;
if the absolute value of θ i 'is less than the first threshold, let θ i' =0; if the absolute value of the theta i ' is not smaller than the first threshold value, enabling the theta i+1＝θi+θi ' to be used as the latest current inclination angle, and correcting the inclination angle theta i ' of the image;
If the absolute value of η i 'is less than the second threshold, let η i' =0; if the absolute value of η i ' is not smaller than the second threshold, η i+1＝ηi+ηi ' is set as the latest current tilt angle, and the rotation angle η i ' is corrected for the image.
The first threshold is δ= 0.00087rad and the second threshold is δ1=0.000087 rad.
In the step (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.
The beneficial effects are that: according to the perspective model theory, the angle error is calculated by utilizing the change of the projection size of the object caused by the rotation of the detector in different directions, so that the accuracy is ensured, and the robustness of angle solving is improved. 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 errors in the calibration process until the calibration precision requirement is met.
Drawings
FIG. 1 is a schematic view of a projection 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 model projection of a flat panel detector with a tilt angle phi;
FIG. 4 is a schematic view of a model projection of a flat panel detector with an inclination angle θ;
FIG. 5 is a schematic view of a flat panel detector with simultaneous tilt angle phi, theta and time scale model projection;
FIG. 6 is a schematic view of a model projection of a flat panel detector with an inclination angle eta.
Detailed Description
The invention is further elucidated below in connection with the drawings and the specific 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 parallel in the calibration die body, the number of each layer of steel balls is 12, the diameter of each steel ball is 2mm, each layer of 12 steel balls are arranged on a standard circular track, the standard circular tracks are 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 is in unit mm; the vertical distance between the upper and lower standard circular tracks is D, in mm. The invention adopts the X-ray source of the cone beam CT to scan the calibration die body and project the calibration die body to the flat panel detector, the X-ray source passes through the calibration die body, the central line of the ray 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 the right-position and side-position shooting images, 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, a sequence image is obtained, as shown in FIG. 1, an X-ray beam emitted by an X-ray source S passes through the calibration phantom, and a projection image is generated on a flat panel detector by projection. The projection of a standard circle formed by the upper steel ball layer and the lower steel ball layer on the plane detector is an ellipse, and the top points of the upper ellipse and the lower ellipse which are farthest from each other are V 0、V1、V2、V3, namely the long axis top points of the ellipses are V 0、V1、V2、V3 respectively.
Fig. 2 is a flowchart of an iterative correction method of an oblique image according to the present invention, as shown in fig. 2, the iterative correction method of an oblique image according to the present invention includes the following steps:
(1) Performing rotary scanning on the calibration phantom through cone beam CT to obtain a projection image of the calibration phantom on the flat panel detector;
(2) Obtaining the center coordinates of each steel ball: carrying out Gaussian filtering and threshold segmentation on a single projection image, determining a segmentation threshold value according to the number of pixel points projected by steel balls, initially extracting the position coordinates of the steel balls through a self-adaptive threshold segmentation algorithm, carrying out edge extraction segmentation on the steel balls according to the position coordinates, obtaining the accurate contour of each steel ball, and carrying out ellipse fitting by adopting a least square method to obtain the center coordinates of each steel ball;
(3) Sequencing the upper steel balls and the lower steel balls obtained in the step (2) according to anticlockwise, and carrying out central symmetry pairing on the upper steel balls and the lower steel balls, wherein the intersection point between the connecting lines of the central symmetry steel balls of each pair is the center of the calibration phantom, and solving the coordinate of the central point of the calibration phantom projected on the flat panel detector according to the intersection point;
(4) Respectively carrying out least square ellipse fitting on the upper and lower steel ball sets to respectively obtain ellipse equations of the upper and lower steel balls, and solving to obtain major axis vertexes and ellipse center coordinates 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 the projection of the calibration phantom onto the flat panel detector when the detector has an inclination angle phi, as shown in FIG. 3, the dashed rectangle indicates the ideal position of the flat panel detector, and the solid rectangle indicates the actual position of the flat panel detector; when the flat panel detector is at an ideal position, the projection of the ray of the X-ray source S passing through the center of the calibration phantom on the flat panel detector is O; an ideal coordinate system is established by taking the point O as an origin, taking an axis which passes through the point O and is parallel to the vertical side of the flat panel detector as a y-axis, taking an axis x-axis which passes through the origin O and is parallel to the horizontal side of the flat panel detector, and taking an axis which passes through the point O and is perpendicular to the x-axis and the y-axis at the same time as a z-axis;
when the inclination angle of the ideal position of the flat panel detector relative to the flat panel detector around the y axis is phi, an included angle exists between an actual coordinate system where the flat panel detector is positioned and the ideal coordinate system, and the actual coordinate system consists of an x φ axis, a y φ axis and a z φ axis, wherein the x φ axis is obtained by rotating the x axis around the y axis by an angle phi in the corresponding direction, the y φ axis is consistent with the y axis, and the z φ axis is obtained by rotating the z axis around the y axis by an angle phi in the corresponding direction;
FIG. 4 is a schematic view of a projection of a calibration phantom onto a flat panel detector when the flat panel detector has an inclination angle θ, as shown in FIG. 4, the dashed rectangle indicates the ideal position of the flat panel detector in an ideal case, and the solid rectangle indicates the actual position of the flat panel detector; when the inclination angle theta of the flat panel detector relative to the ideal position of the flat panel detector around the x-axis is equal to the inclination angle theta of the flat panel detector, an included angle exists between an actual coordinate system where the flat panel detector is actually located and an ideal coordinate system where the ideal position of the flat panel detector is located, the actual coordinate system consists of an x θ axis, a y θ axis and a z θ axis, wherein the x θ axis is consistent with the x axis, the y θ axis is obtained by rotating the y axis around the x-axis, and the z θ axis is obtained by rotating the z axis around the x-axis;
FIG. 5 is a schematic diagram of a model projection of a detector in which the inclination angles phi and theta exist simultaneously, and the current inclination angle values phi i and theta i of the inclination angles phi and theta of the flat panel detector are calculated according to formulas (1) and (2) respectively;
Wherein V 0φ、V1φ is the vertex of the long axis of an ellipse formed by the projection of a circle formed by the upper layer of steel balls on an actual flat panel detector under the condition that the inclination angle of the ideal position of the flat panel detector relative to the flat panel detector around the y axis is phi i; v 2φ、V3φ is the vertex of the long axis of an ellipse formed by the projection of a circle formed by the lower layer steel balls on an actual flat panel detector under the condition that the inclination angle of the flat panel detector around the y axis relative to the ideal position of the flat panel detector is phi i; v 0θ、V1θ is the vertex of the long axis of an ellipse formed by the projection of a circle formed by the upper layer of steel balls on an actual flat panel detector under the condition that the inclination angle of the flat panel detector around the x axis relative to the ideal position of the flat panel detector is theta i; v 2θ、V3θ is the vertex of the long axis of an ellipse formed by the projection of a circle formed by the lower layer steel balls on an actual flat panel detector under the condition that the inclination angle of the flat panel detector around the x axis relative to the ideal position of the flat panel detector is theta i; d 1 is the distance from the X-ray source S to the center of the calibration phantom, d 2 is the distance from the X-ray source S to the flat panel detector, d 1、d2 is a design value during the primary calculation, and d 1、d2 is the actual distance calculated in the formula (7) during the later iteration;
The projection image is subjected to inclination correction by the calculated existing inclination angles phi i and theta i of the current flat panel detector, ellipse parameters in the corrected image are obtained, the current value eta i of the inclination 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 parameters, and a model projection schematic diagram with the existing inclination angle eta of the flat panel detector is shown in fig. 6, wherein the model projection schematic diagram is specifically calculated as follows:
As shown in fig. 3, the intersection point of the straight lines where V 0φV1φ and V 2φV3φ are located is P φ,Pφ, and the perpendicular point on the y φ axis is P a, which can be calculated by formula (3):
Wherein d Oa is the distance from the point P a to the origin O, d O1 is the distance from the center of the upper ellipse to the origin O, d O2 is the distance from the center of the lower ellipse to the origin O, a 1、b1 is the coefficient of the upper ellipse equation, and a 2、b2 is the coefficient of the lower ellipse equation;
Calculating a current value eta i of an inclination eta of the flat panel detector around the z axis relative to an ideal position of the flat panel detector according to the formula (4);
η＝[d1aA(x,L1)+d2aA(x,L2)]/d12 (4)
Wherein L 1、L2 represents straight lines where V 0φV1φ and V 2φV3φ are located, respectively, A (x, L 1) and A (x, L 2) are respectively And/>The included angle between the straight line and the x axis; d 1a、d2a represents the distance from the center of the upper and lower ellipse to the point P a, and d 12 represents the distance between the centers of the upper and lower ellipse;
(6) Calculating an actual distance d 1 from the X-ray source S to the center of the calibration phantom and an actual distance d 2 from the X-ray source S to the flat panel detector after the flat panel detector is corrected;
Calculating a rotation center Q (namely the center of the calibration phantom) and coordinate values of the X-ray source S in a world coordinate system and an actual coordinate system where the flat panel detector is actually positioned according to current inclination angle values phi i、θi and eta i of the flat panel detector relative to the ideal position of the flat panel detector, and calculating an actual distance d 1 from the X-ray source S to the center of the calibration phantom and an actual distance d 2 from the X-ray source S to the flat panel detector according to the coordinate values; the method comprises the following steps:
(61) Calculating coordinate values of the X-ray source S in an actual coordinate system where the current flat panel detector is located according to the formula (5):
meanwhile, calculating to obtain the coordinate value of the rotation center Q in the actual coordinate system where the current flat panel detector is located by the aid of the method (6):
(62) The actual distance d 1 of the X-ray source S to the center of the calibration phantom and the actual distance d 2 of the X-ray source S to the flat panel detector are calculated according to equation (7):
(7) Calculating a corrected inclination angle phi i′、θi、ηi 'between the flat panel detector corrected by the current value phi i、θi、ηi and the ideal position of the flat panel detector, and comparing the inclination angle phi i′、θi、ηi' with a set threshold value;
Calculating the new actual distance d 1 from the X-ray source S to the center of the calibration phantom and the actual distance d 2 from the X-ray source S to the flat panel detector according to the formula (7), and respectively calculating the corrected inclination angle phi i′、θi' of the corrected flat panel detector according to the formulas (1) and (2);
Setting a first threshold δ= 0.00087rad, and if the absolute value of Φ i 'is smaller than the first threshold, making Φ i' =0; if the absolute value of phi i ' is not less than the first threshold value, phi i+1＝φi+φi ' is made, and the inclination angle phi i ' of the image is corrected;
If the absolute value of θ i 'is less than the first threshold, let θ i' =0; if the absolute value of the theta i ' is not smaller than the first threshold value, enabling the theta i+1＝θi+θi ' to be adopted, and correcting the image inclination angle theta i ';
Calculating through the formula (4) to obtain the inclination angle eta i ' of the corrected flat panel detector, setting a second threshold delta 1 =0.000087 rad, and if the absolute value of eta i ' is smaller than the second threshold, setting eta i ' =0; if the absolute value of eta i ' is not less than the second threshold, enabling eta i+1＝ηi+ηi ' and correcting the image inclination angle eta i ';
(8) Judging whether phi i′、θi′、ηi' is 0 or not, if so, ending the cycle, and calculating a calibration parameter by taking the current inclination angle as an included angle between the actual flat panel detector and the ideal detector position and outputting the calibration parameter, wherein the calibration parameter is a rotation center calculated according to the current inclination angle, a coordinate value of a ray source in a flat panel detector coordinate system, an actual distance d 1 between the ray source and the center of a calibration die body and an actual distance d 2 between the ray source and the flat panel detector; otherwise, return to step (6) with phi i+1、θi+1、ηi+1 as the latest current value phi i、θi、ηi.
In other embodiments, the number of steel balls in each layer can be set to be 4-16, the diameters 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 steel balls have equal diameters.
In the present invention, specific values of the first threshold value and the second threshold value may be set according to actual conditions.
According to the perspective model theory, the angle error is calculated by utilizing the change of the projection size of the object caused by the rotation of the flat panel detector in different directions, so that the accuracy is ensured and the robustness of angle solving is improved. Meanwhile, in order to eliminate the mutual image problem among different parameters, the invention adopts an iterative calibration and successive correction mode to gradually correct the geometric parameter errors until the calibration precision requirement is met.
The preferred embodiments of the present invention have been described in detail above, but the present invention is not limited to the specific details of the above embodiments, and various equivalent changes (such as number, shape, position, etc.) may be made to the technical solution of the present invention within the scope of the technical concept of the present invention, and these equivalent changes all fall within the scope of the present invention.
Claims (10)
1. An iterative correction method for an oblique image is characterized in that: the method comprises the following steps:
(1) Scanning the calibration die body through CT to obtain a projection image of the calibration die body on the flat panel detector, and obtaining the center coordinates of steel balls in the projection image, wherein the steel balls are arranged in two layers and are arranged in parallel in the calibration die body, and each layer of steel balls form a standard circular track;
(2) Respectively calculating equations of ellipses formed on the projection images, and solving corresponding ellipse parameters according to the equations;
(3) Calculating a current inclination angle of the current flat panel detector relative to an ideal position of the flat panel detector according to the design parameters of the calibration die body, the ellipse parameters obtained in the step (2) and the distance between the radiation source and the flat panel detector, wherein the current inclination angle is specifically as follows:
The flat panel detector has a tilt angle phi about an axis parallel to its vertical side and a tilt angle theta about an axis parallel to its horizontal side with respect to its ideal position;
Respectively calculating corresponding dip angles phi i and theta i according to formulas (1) and (2);
wherein R is the radius of a circular track formed by each layer of steel balls, and D is the vertical distance between circular tracks formed by two layers of steel balls which are arranged in parallel; v 0φ、V1φ is the vertex of the long axis of the ellipse formed by the projection of the circle formed by the upper steel ball on the actual flat panel detector under the condition that the inclination angle of the flat panel detector relative to the ideal position of the flat panel detector around the axis parallel to the vertical edge of the flat panel detector is phi i; v 2φ、V3φ is the vertex of the long axis of the ellipse formed by the projection of the circle formed by the lower steel ball on the actual detector under the condition that the inclination angle of the flat panel detector relative to the ideal position of the flat panel detector around the axis parallel to the vertical edge of the flat panel detector is phi i; v 0θ、V1θ is the vertex of the long axis of the ellipse formed by the projection of the circle formed by the upper layer steel ball on the actual detector under the condition that the inclination angle of the flat panel detector relative to the ideal position of the flat panel detector around the axis parallel to the horizontal edge of the flat panel detector is theta i; v 2θ、V3θ is the vertex of the long axis of the ellipse formed by the projection of the circle formed by the lower steel ball on the actual detector under the condition that the inclination angle of the flat panel detector relative to the ideal position of the flat panel detector around the axis parallel to the horizontal edge of the flat panel detector is theta i; d 2 is the distance from the radiation source to the flat panel detector, and the initial value is a design value;
Performing inclination correction on the projection image by using corresponding inclination angles phi i and theta i to obtain ellipse parameters in the corrected image, and calculating an inclination angle eta i of the flat panel detector diffraction line source taking the connection line of the flat panel detector projection and the ray source as an axis according to the ellipse parameters;
The calculation of the inclination angle eta i specifically comprises the following steps:
The intersection point of the straight lines where V 0φV1φ and V 2φV3φ are located is P φ,Pφ, the vertical point of the flat panel detector is P a, and the distance d Oa from the point P a to the origin O is calculated by the formula (3):
Wherein d O1 is the distance from the center of the upper ellipse to the origin O, d O2 is the distance from the center of the lower ellipse to the origin O, a 1、b1 is the coefficient of the upper ellipse equation, and a 2、b2 is the coefficient of the lower ellipse equation;
Calculating according to the formula (4), obtaining the inclination angle eta i of the flat panel detector diffraction line source taking the projection of the flat panel detector and the line connecting with the ray source as the axis:
η＝[d1aA(x,L1)+d2aA(x,L2)]/d12 (4)
Wherein L 1、L2 represents straight lines in which V 0φV1φ and V 2φV3φ are located, respectively, and A (x, L 1) and A (x, L 2) are respectively And/>The included angle between the straight line and the axis parallel to the horizontal side of the flat panel detector; d 1a、d2a represents the distance from the center of the upper and lower ellipse to the point P a, and d 12 represents the distance between the centers of the upper and lower ellipse;
(4) Correcting the inclination angle of the image by using the current inclination angle of the flat panel detector, and calculating corrected ellipse parameters and the actual distance from the ray source to the flat panel detector;
(5) Obtaining ellipse parameters and actual distance between the ray source and the flat panel detector according to the design parameters of the calibration die body and the step (4), and calculating to obtain corrected inclination angles of the flat panel detector relative to the ideal positions of the flat panel detector after correction;
(6) Judging whether the inclination angle corrected in the step (5) is smaller than a preset threshold value, if so, ending the cycle; otherwise, the sum of the corrected inclination angle obtained in the step (5) and the current inclination angle in the step (4) is used as the latest current inclination angle, and the step (4) is returned.
2. The iterative correction method of oblique images according to claim 1, characterized in that: in the step (1), the number of steel balls in each layer in the calibration die body is set to be 4-16, and the diameters of the steel balls are set to be 1-5 mm.
3. The iterative correction method of oblique images according to claim 2, characterized in that: in the calibration die body, the diameter of each steel ball is the same.
4. A method of iterative correction of oblique images as claimed in claim 3, characterized in that: the number of steel balls in each layer is set to 12, and the diameter of each steel ball is set to 2mm.
5. The iterative correction method of oblique images 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 of oblique images according to claim 1, characterized in that: in the step (1), the center coordinates of the steel balls are obtained specifically as follows:
And carrying out Gaussian filtering and threshold segmentation on the image, determining a segmentation threshold value according to the number of pixel points projected by the steel balls, initially extracting the position coordinates of the steel balls through a self-adaptive threshold segmentation algorithm, carrying out edge extraction segmentation on the steel balls according to the position coordinates, obtaining the accurate contour of each steel ball, and carrying out ellipse fitting by adopting a least square method to obtain the center coordinates of each steel ball.
7. The iterative correction method of oblique images according to claim 1, characterized in that: the step (5) is to calculate the corrected ellipse corresponding parameter and the distance between the ray source and the actual flat panel detector specifically as follows:
calculating the coordinate value of the ray source S in the actual coordinate system where the current flat panel detector is located according to the formula (5):
The actual distance d 2 from the source S to the actual flat panel detector is calculated according to the following equation:
8. The iterative correction method of oblique images according to claim 7, characterized in that: the step (6) specifically comprises the following steps:
calculating the actual distance d 2 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 i'、θi'、ηi' according to the ellipse parameters obtained in the step (4) and the combination formulas (1), (2) (3) and (4);
If the absolute value of phi i 'is less than the first threshold, let phi i' =0; if the absolute value of phi i ' is not less than the first threshold value, phi i+1＝φi+φi ' is taken as the latest current inclination angle, and the inclination angle phi i ' of the image is corrected;
If the absolute value of θ i 'is less than the first threshold, let θ i' =0; if the absolute value of the theta i ' is not smaller than the first threshold value, enabling the theta i+1＝θi+θi ' to be used as the latest current inclination angle, and correcting the inclination angle theta i ' of the image;
If the absolute value of η i 'is less than the second threshold, let η i' =0; if the absolute value of η i ' is not smaller than the second threshold, η i+1＝ηi+ηi ' is set as the latest current tilt angle, and the rotation angle η i ' is corrected for the image.
9. The iterative correction method of oblique images according to claim 8, characterized in that: the first threshold is set to δ= 0.00087rad and the second threshold is set to δ 1 =0.000087 rad.
10. The iterative correction method of an oblique image according to any one of claims 1 to 9, characterized in that: in the step (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 CN113876346A (en) | 2022-01-04 |
| CN113876346B true 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) |
Families Citing this family (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115105757B (en) * | 2022-05-24 | 2025-07-18 | 沈阳东软智睿放疗技术有限公司 | Medical die body, and equipment precision detection method and system based on medical die body |
| 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 |
| 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 |
Citations (2)
| 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 |
| CN108122203A (en) * | 2016-11-29 | 2018-06-05 | 上海东软医疗科技有限公司 | A kind of bearing calibration of geometric parameter, device, equipment and system |
Family Cites Families (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| DE602007010192D1 (en) * | 2006-01-26 | 2010-12-16 | Toshiba Kk | X-ray CT device, as well as methods for phantom assignment |
| EP2884897A1 (en) * | 2012-08-20 | 2015-06-24 | orangedental GmbH & Co. KG | Geometric characterization and calibration of a cone-beam computer tomography apparatus |
| GB2520711B (en) * | 2013-11-28 | 2018-06-20 | Nikon Metrology Nv | Calibration apparatus and method for computed tomography |
-
2021
- 2021-11-16 CN CN202111352681.4A patent/CN113876346B/en active Active
Patent Citations (2)
| 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 |
| 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 机设计及其关键技术研究;张震;工程科技Ⅱ辑(第06期);1-71 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN113876346A (en) | 2022-01-04 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN113876346B (en) | Iterative correction method for inclined image | |
| CN109598762B (en) | High-precision binocular camera calibration method | |
| WO2018218611A1 (en) | Geometric parameter determination method for cone beam computed tomography system | |
| CN103735282B (en) | A kind of cone-beam CT system detector geometric correction device and bearing calibration thereof | |
| CN106667512B (en) | Geometric correction method of X-ray imaging device and breast tomography device | |
| CN112603346B (en) | A Detector Deflection Correction Method Based on Marker Imaging | |
| CN112762910B (en) | Short-measuring-range correction calibration method suitable for laser scanner | |
| CN110084855B (en) | Improved CBCT geometric parameter calibration method | |
| CN117392211B (en) | BGA element rapid identification positioning method and system and storage medium | |
| CN114170284A (en) | Multi-view point cloud registration method based on active landmark projection assistance | |
| CN106706675A (en) | Correction method based on computed laminography (CL) system | |
| CN110490941B (en) | Telecentric lens external parameter calibration method based on normal vector | |
| CN112971984B (en) | Coordinate registration method based on integrated surgical robot | |
| CN110645928A (en) | Space coordinate positioning method of three-dimensional scanner | |
| CN115100289A (en) | Camera sensor installation state calibration method | |
| CN118014911A (en) | Correction phantom, CBCT geometric correction method, system and device based on offset detector | |
| CN112132903B (en) | Coordinate system calibration method and system for visual system and multi-axis motion system | |
| CN203763103U (en) | Geometric correction device of detector of cone beam CT (computed Tomography) system | |
| CN117571719A (en) | A method for locating vehicle body paint defects based on pixel back projection | |
| CN110501360B (en) | Standard device for correcting pose of micro CT system and implementation method | |
| CN115018921B (en) | Camera calibration plate and calibration method thereof | |
| CN117974826B (en) | Cone beam CT geometric parameter calibration method, device, equipment and storage medium | |
| CN119934976A (en) | Calibration method of microstructure light three-dimensional measurement system based on standard spherical surface | |
| CN117853593B (en) | Linear array camera calibration method based on two-dimensional code | |
| CN115082543B (en) | Laser correction method |
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 | ||
| 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. |
|
| GR01 | Patent grant |