CN111968164A - Automatic implant registration and positioning method based on biplane X-ray tracking - Google Patents
Automatic implant registration and positioning method based on biplane X-ray tracking Download PDFInfo
- Publication number
- CN111968164A CN111968164A CN202010836044.3A CN202010836044A CN111968164A CN 111968164 A CN111968164 A CN 111968164A CN 202010836044 A CN202010836044 A CN 202010836044A CN 111968164 A CN111968164 A CN 111968164A
- Authority
- CN
- China
- Prior art keywords
- ray
- obj
- image
- coordinate system
- implant
- 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.)
- Pending
Links
- 239000007943 implant Substances 0.000 title claims abstract description 118
- 238000000034 method Methods 0.000 title claims abstract description 50
- 238000002922 simulated annealing Methods 0.000 claims abstract description 4
- 239000011159 matrix material Substances 0.000 claims description 91
- 230000002285 radioactive effect Effects 0.000 claims description 17
- 230000009466 transformation Effects 0.000 claims description 14
- 238000011156 evaluation Methods 0.000 claims description 10
- 230000033001 locomotion Effects 0.000 claims description 10
- 230000005855 radiation Effects 0.000 claims description 8
- 238000002513 implantation Methods 0.000 claims description 7
- 241001394244 Planea Species 0.000 claims description 6
- 238000000605 extraction Methods 0.000 claims description 6
- 238000000137 annealing Methods 0.000 claims description 4
- 230000009467 reduction Effects 0.000 claims description 4
- 238000010521 absorption reaction Methods 0.000 claims description 3
- 230000010339 dilation Effects 0.000 claims description 3
- 230000008961 swelling Effects 0.000 claims description 3
- 238000004422 calculation algorithm Methods 0.000 abstract description 4
- 238000004458 analytical method Methods 0.000 abstract description 3
- 210000000629 knee joint Anatomy 0.000 description 10
- 210000000988 bone and bone Anatomy 0.000 description 7
- 239000002184 metal Substances 0.000 description 7
- 238000002591 computed tomography Methods 0.000 description 6
- 238000011883 total knee arthroplasty Methods 0.000 description 6
- 238000006073 displacement reaction Methods 0.000 description 5
- 238000005516 engineering process Methods 0.000 description 5
- 230000008569 process Effects 0.000 description 5
- 238000010586 diagram Methods 0.000 description 4
- 230000000694 effects Effects 0.000 description 4
- 230000002980 postoperative effect Effects 0.000 description 4
- 238000013527 convolutional neural network Methods 0.000 description 3
- 238000003384 imaging method Methods 0.000 description 3
- 210000003127 knee Anatomy 0.000 description 3
- 210000002303 tibia Anatomy 0.000 description 3
- 210000000689 upper leg Anatomy 0.000 description 3
- NIXOWILDQLNWCW-UHFFFAOYSA-N acrylic acid group Chemical group C(C=C)(=O)O NIXOWILDQLNWCW-UHFFFAOYSA-N 0.000 description 2
- 210000003484 anatomy Anatomy 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 2
- 230000005540 biological transmission Effects 0.000 description 2
- 230000007547 defect Effects 0.000 description 2
- 238000003708 edge detection Methods 0.000 description 2
- 238000001727 in vivo Methods 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 201000008482 osteoarthritis Diseases 0.000 description 2
- 208000003947 Knee Osteoarthritis Diseases 0.000 description 1
- 239000004698 Polyethylene Substances 0.000 description 1
- 238000004364 calculation method Methods 0.000 description 1
- 239000000919 ceramic Substances 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 230000002950 deficient Effects 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 239000000284 extract Substances 0.000 description 1
- 238000002594 fluoroscopy Methods 0.000 description 1
- 230000002068 genetic effect Effects 0.000 description 1
- 230000036541 health Effects 0.000 description 1
- 208000015181 infectious disease Diseases 0.000 description 1
- 210000003041 ligament Anatomy 0.000 description 1
- 239000000463 material Substances 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- 230000000399 orthopedic effect Effects 0.000 description 1
- -1 polyethylene Polymers 0.000 description 1
- 229920000573 polyethylene Polymers 0.000 description 1
- 238000012805 post-processing Methods 0.000 description 1
- 230000000750 progressive effect Effects 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 230000003068 static effect Effects 0.000 description 1
- 238000001356 surgical procedure Methods 0.000 description 1
- 238000003325 tomography Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/30—Determination of transform parameters for the alignment of images, i.e. image registration
- G06T7/33—Determination of transform parameters for the alignment of images, i.e. image registration using feature-based methods
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/70—Determining position or orientation of objects or cameras
- G06T7/73—Determining position or orientation of objects or cameras using feature-based methods
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/10—Image acquisition modality
- G06T2207/10116—X-ray image
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/30—Subject of image; Context of image processing
- G06T2207/30004—Biomedical image processing
- G06T2207/30008—Bone
Landscapes
- Engineering & Computer Science (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Apparatus For Radiation Diagnosis (AREA)
Abstract
The invention provides an implant automatic registration positioning method based on biplane X-ray tracking, which comprises the steps of constructing a virtual biplane X-ray system, extracting and weighting expansion of an implant edge in an X-ray image, rapidly projecting and extracting an implant geometric model, and searching an optimal matching position. Determining the relative position of the image planes of the two X-ray machines according to a locator, thereby establishing a virtual double-plane X-ray system in a computer system; the weighted expansion edge extracted from the X-ray images shot by the two X-ray machines and the projected edge of the implant model are utilized to calculate the matching degree of the implant model at any spatial position, and the spatial position of the implant at the highest matching degree can be rapidly obtained through a simulated annealing algorithm. The method can realize automatic tracking of the implant in the body space position based on the moving images continuously obtained by two X-ray machines, and the method is used as a reference to carry out the kinematics research and the joint contact analysis of the implant, thereby greatly reducing the manpower and time expenditure of manual registration.
Description
Technical Field
The invention relates to the fields of in-vivo kinematic tracking, X-ray image post-processing and 2D-3D image registration of metal implants, in particular to an automatic implant registration and positioning method based on biplane X-ray tracking.
Background
Total Knee Arthroplasty (TKA) is a procedure used in clinical orthopedics to treat end-stage osteoarthritis in the Knee. TKA is mainly to remove defective bone tissues at the distal femur and the proximal tibia, and to reconstruct a joint using an artificial implant made of metal, ceramic, high polyethylene, etc., thereby eliminating pain in the knee joint of a patient. With the development of the design principle, materials, configuration and the like of the implant, the TKA postoperative patient can recover the motion function to a certain extent, and the quality of life and the independent living capacity are improved. The preoperative doctor can comprehensively evaluate factors such as the progressive degree of knee osteoarthritis, the function and integrity of cruciate ligaments and the age of a patient, so as to select a surgical formula and a knee implant which are most suitable for the patient and determine the configuration of the implant to be used. Clinical reports indicate that 15-30% of patients are not satisfied with the effect after TKA [1], the proportion of patients with the effective service life of the TKA postoperative implant being 10-15 years is more than 90%, but the rest patients needing implant revision bring 270 billion dollars of burden to the health system of the United states [2 ]. Factors that mainly cause the need for secondary revision surgery for knee implants include infection, implant loosening, sex, age, etc. [1-3 ]. The implant implantation position is determined by a doctor in an operation according to the joint state in a non-load bearing state, a postoperative patient needs to carry out a large number of rehabilitation daily life actions such as sitting down, getting up, walking, going up and down stairs and the like, the normal anatomical structure of the knee joint cannot be completely recovered due to the fact that the spatial motion and the joint surface contact position of the knee joint implant are different from the static implantation position in the operation, and therefore the implant can be loosened and fails, and secondary operations are needed. Therefore, the evaluation of the in-vivo dynamic function of the knee joint is beneficial to determining the optimal implantation position of the knee joint implant, improving the configuration of the implant and making a personalized postoperative rehabilitation plan.
The biplane dynamic X-ray tracking technology (DFIS) is a technology for synchronously recording joint motion images by using two orthogonal X-ray machines and accurately measuring the spatial position and the motion function of a joint [4-6 ]. The method comprises the steps of synchronously recording kinematic image data of joints by using two mobile X-ray machines, acquiring data by using a Computed Tomography (CT) technology and a Magnetic Resonance (MR) technology, establishing a 3D model of the personalized bone and knee joint implant, and matching the 3D bone model with a 2D medical image in a virtual biplane system based on a point light source projection principle, so as to obtain the spatial relative position of the joints for subsequent kinematic analysis and other parameter analysis. The precision of the biplane dynamic X-ray tracking technology for measuring the displacement and the rotation angle of the dynamic movement of the knee joint is 0.2mm and 0.2 degrees respectively [5 ]. However, in the current virtual biplane system, an operator is required to manually adjust the spatial position of the 3D model to realize the 2D-3D registration process, the manual registration process is subjective, the matching results of different operators may be different, and a lot of time is consumed. A skilled operator can quickly adjust the position of the 3D model to the vicinity of the exact position so that the projected contour and the 2D image contour approximately match (takes 10-20 seconds), but to obtain the exact 3D model position it takes 3-5 minutes to fine-tune the single model position, whereas the knee joint is made up of femur and tibia, the model registration of a set of 2D images takes more than 5 minutes for the operator, and it takes 300-400 minutes to fully register a set of longer-cycle motion data (30-40 sets of 2D images). Therefore, the method develops an implant automatic registration positioning program based on biplane X-ray tracking, improves the repeatability and objectivity of the registration result, greatly reduces the manpower and time expenditure in the registration process, and is very valuable.
At present, there are 2D-3D automatic registration and positioning methods based on single plane X-rays, such as mahfuz and other automatic registration methods based on contour feature matching for 2D single plane X-rays and 3D metal implant models [7], Tsai and other automatic registration methods based on digital Reconstructed radiographic images (DRRs) for 2D single plane X-rays and 3D digital tomographic imaging data (CT) [8], Miao and other real-time automatic registration methods based on convolutional neural network for 2D single plane X-rays and 3D metal implant models [9 ]. The basic principle of the automatic registration methods is that in a virtually-constructed single-plane projection system, a virtual 2D plane projection image is generated according to the relative positions of a 3D target object, a radioactive source and a projection plane, the similarity degree of the virtual 2D plane projection image and a 2D single-plane X-ray image recorded in actual transmission is calculated, and the spatial position of the 3D target object is adjusted by using an optimization algorithm so that the similarity degree of the projected 2D image and the 2D single-plane X-ray image recorded in actual transmission is the highest. There are 3 key technical issues in technical implementation:
(1) method of generating a virtual 2D planar projection image.
Firstly, a 3D skeleton and metal implantation model is a 3D triangular patch set which is formed by a series of nodes at specific positions and the interconnection relation between the nodes in a computer system, the nodes and the triangular patches can be directly projected onto a projection plane according to the basic theorem of perspective projection, and a 2D binarization mask (mask) is generated to represent the projection result of the 3D model [7, 9 ].
CT is a 3D volume data set composed of a series of volume data units with the same space size, the anatomical structures and relative positions of bones and metal implants under no load can be reconstructed, and a DRRs method [10] can be used for generating a 2D virtual X-ray perspective projection image [8] based on the relation among a virtual radiation source, CT data and a projection plane.
(2) And (5) a similarity evaluation method.
[7] A similarity evaluation method based on gray level distribution and contour features. The virtual 2D plane projection image obtained based on the 3D model projection is a binary image (with projection as 1, without projection as 0), the actually recorded 2D single plane X-ray image is an 8-bit grayscale image (grayscale range 0-255), and the normalized grayscale matching degree can be calculated using the following formula:
wherein G (X, y) is a gradient matrix of the actually recorded 2D single-plane X-ray image, and H (X, y) is a binary matrix of the virtual 2D plane projection image profile. The overall matching degree can be calculated by the weighted summation of the gray matching degree and the contour matching degree.
Extracting the geometric contour of the virtual 2D plane projection image, expressing the geometric contour by a binary matrix (the contour is marked as 1, and the non-contour is marked as 0), calculating the gradient value of the actually recorded 2D single plane X-ray image, wherein the gradient value of the contour is high, the gradient value of the non-contour is low, and the normalized contour matching degree can be calculated by using the following formula:
wherein J (X, y) is a gray matrix of the actually recorded 2D single-plane X-ray image, and J (X, y) is a binary matrix of the virtual 2D plane projection image.
② an evaluation method based on gradient Matching degree and Weighted contour Matching degree (WEMS) [8 ]. The virtual 2D plane projection image projected by the DRR method and the actually recorded 2D single plane X-ray image are both gray level images, the gradient values can be respectively calculated, and the gradient value matching degree is calculated by using the following formula:
wherein If1(X, y) is the gray matrix of the actual recorded 2D single-plane X-ray image, IDRR(x, y) is the DRR-generated virtual 2D plane projection image gray-scale matrix.
And (3) extracting the edge of the actually recorded 2D single-plane X-ray image by using a Canny edge detection operator, generating a binary image (the edge is recorded as 1, and the non-edge is recorded as 0), and performing weighted expansion on the binary edge image by using a binary image expansion operation. Similarly, extracting the edge of the virtual 2D plane projection image by using a Canny edge detection operator, generating a binary image (the edge is recorded as 1, and the non-edge is recorded as 0), and performing weighted assignment on the binary edge image according to the length of the edge, wherein the longer the edge is, the higher the weight is. The weighted edge match metric (WEMS) is calculated using the following formula:
Ereg(x，y)＝LDRR(x，y)Bf1(x，y).
wherein B isf1(X, y) is the expanded, actually recorded 2D monoplane X-ray image edge weight matrix, LDRRAnd (x, y) is an edge weight matrix of the virtual 2D plane projection image after being assigned according to the length.
(3) And optimizing the calculation method.
Simulation Annealing algorithm (SA) [7 ];
genetic Algorithm (GA) [8 ];
(iii) Convolutional Neural Network (CNN) [9 ].
The above prior art has the following disadvantages:
(1) the conventional biplane dynamic X-ray imaging system requires manual registration, the manual registration process is subjective, the matching result may be different for different operators, and in addition, a great amount of time is consumed. A skilled operator can quickly adjust the position of the 3D model to the vicinity of the exact position so that the projected contour and the 2D image contour approximately match (takes 10-20 seconds), but to obtain the exact 3D model position it takes 3-5 minutes to fine-tune the single model position, whereas the knee joint is made up of femur and tibia, the model registration of a set of 2D images takes more than 5 minutes for the operator, and it takes 300-400 minutes to fully register a set of longer-cycle motion data (30-40 sets of 2D images).
(2) The 2D-3D automatic registration positioning method based on the single plane X-ray has defects. Essentially, single plane X-rays are a single view imaging system. By adjusting the position of the 3D object in the projection plane direction, the virtual 2D projection image can be translated significantly, which significantly affects the matching degree of the virtual 2D image and the actual 2D image, but by adjusting the position of the 3D object in the out-of-plane direction, the virtual 2D projection image can be enlarged or reduced in size, and the enlargement or reduction ratio is affected by the out-of-plane displacement of the 3D object and the distance between the radiation source and the projection plane. When the radioactive source is close to the projection plane, the small displacement of the 3D target object in the out-of-plane direction can cause larger size change of the virtual 2D projection image, but the space between the radioactive source and the projection plane is smaller, and the perspective space is greatly limited; conversely, when the radiation source is far from the projection plane, the perspective space is guaranteed, but small displacements in the "out-of-plane" direction of the 3D object have little effect on the virtual 2D projection image. In order to ensure that the knee joint can be imaged by X-ray fluoroscopy as much as possible when the function of the subject moves, the distance between a radioactive source of an X-ray machine and a projection plane is kept about 1 meter, which results in that a similarity evaluation function in the 2D-3D registration process is insensitive to the displacement of a 3D target object in the out-of-plane direction, so that the 2D-3D registration based on a single-plane image has a larger error in the out-of-plane direction. The automatic registration method of the 2D single-plane X-ray and 3D metal implant model based on contour feature matching, such as Mahfouz, has 1.4mm out-of-plane direction error [7 ]; tsai and other automatic registration methods based on 2D single-plane X-ray and 3D digital tomography data (CT) of a digital reconstructed image have a maximum 7.4mm out-of-plane direction error [8] when tibial bone registration is carried out, because when the 2D single-plane X-ray image of the bone extracts the edge, the bone edge of a non-target area is extracted, and the result of a similarity evaluation function is interfered.
The documents referred to above are as follows: the reference is given by number, which is the number corresponding to the preceding of the following documents.
1.Ilana N.Ackerman，Megan A.Bohensky，Ella Zomer，et al.，The projected burden of primary total knee and hip replacement for osteoarthritis in Australia to the year 2030.Bmc Musculoskeletal Disorders，2019.20(1).
2.Douglas A.Dennis，Richard D.Komistek，and Mohamed R.Mahfouz，In Vivo Fluoroscopic Analysis Of Fixed-Bearing Total Knee Replacements.Clinical Orthopaedics&Related Research，2003.410：p.114-130.
3.Francesco Castagnini，Alessandra Sudanese，Barbara Bordini，et al.，Total Knee Replacement in Young Patients：Survival and Causes of Revision in A Registry Population.Journal of Arthroplasty，2017：p.S0883540317304916.
4.Tsung-Yuan Tsai，MingHanLincoln Liow，Guoan Li，et al.，Bi-Cruciate Retaining Total Knee Arthroplasty Does Not Restore Native Tibiofemoral Articular Contact Kinematics During Gait.Journal of Orthopaedic Research，2019.37(9).
5.Jeffrey Bingham and Guoan Li，An Optimized Image Matching Method for Determining In-Vivo TKA Kinematics with a Dual-Orthogonal Fluoroscopic Imaging System.Journal of Biomechanical Engineering，2006.128(4)：p.588.
6.Christophe Sauret，Hélène Pillet，Wafa Skalli，and Morgan Sangeux，On the use of knee functional calibration to determine the medio-lateral axis of the femur in gait analysis：Comparison with EOS biplanar radiographs as reference.Gait&Posture，2016.50.
7.Mohamed R.Mahfouz，William A.Hoff，Richard D.Komistek，and Douglas A.Dennis，A robust method for registration of three-dimensional knee implant models to two-dimensional fluoroscopy images.IEEE Transactions on Medical Imaging，2003.22(12)：p.p.1561-1574.
8.Tsung-Yuan Tsai，Tung-Wu Lu，Chung-Ming Chen，et al.，A volumetric model-based 2D to 3D registration method for measuring kinematics of natural knees with single-plane fluoroscopy.Medical Physics，2010.37(3).
9.Shun Miao，Z.Jane Wang，and Liao Rui，A CNN Regression Approach for Real-Time 2D/3D Registration.IEEE Transactions on Medical Imaging，2016.35(5)：p.1352-1363.
10.L.Lemieux，R.Jagoe，D.R.Fish，et al.，A patient-to-computed-tomography image registration method based on digitally reconstructed radiographs.Medical Physics，1994.21(11)：p.1749-1760.
Disclosure of Invention
Aiming at the defects in the prior art, the invention aims to provide an implant automatic registration positioning method based on biplane X-ray tracking. The technical scheme of the invention is as follows:
an implant automatic registration positioning method based on biplane X-ray tracking comprises the following steps:
s1: constructing a virtual biplane X-ray system, comprising:
s11: preparing two X-ray machines: the X-ray machine comprises a No. 1X-ray machine and a No. 2X-ray machine, and is used for shooting X-ray images of implants in a patient body;
s12: in a computer system, determining the relative positions of two X-ray machines;
s13: placing an implant model in a computer system in a public area projected by two X-ray machines;
s2: extracting the edge of the implant in the X-ray image of the implant, and performing weighted expansion;
s3: performing fast projection and edge extraction on the implant model;
s4: and performing virtual biplane automatic registration based on the edge features obtained in the steps S2 and S3.
Optionally, the S11 further includes:
two X-ray machines were placed in the laboratory space: no. 1X-ray machine, No. 2X-ray machine make the patient who accomplishes the implant implantation carry out the function motion under the operator guide, and two X-ray machines shoot the implant and implant the joint simultaneously, obtain the X image respectively.
Optionally, the step S12 further includes:
s121: a positioner is arranged in a common projection area of two X-ray machines, 4 lead points are arranged in the middle of the positioner and distributed to form a square, the lead points have X-ray absorption and can be represented by low-gray circular points in an X-ray picture;
s122: the number 1X-ray machine is marked as F1, and a rigid orthogonal space coordinate system (R) is established by taking an X-ray image space shot by F1 as a reference1，V1) Wherein R is1Is a 3X 3 size matrix, each column represents a coordinate axis unit vector of the coordinate system of the X-ray image captured by F1; v1The image coordinate system is a 3 multiplied by 1 size column vector, represents the original point position of the F1 image coordinate system, and takes the image center as the original point position of the F1 image coordinate system; when the F1 image space is taken as the system coordinate system, it can be known that:
s123: the No. 2X-ray machine is marked as F2, and a rigid orthogonal space coordinate system (R) is established by taking an X-ray image space shot by F2 as a reference2，V2) Wherein R is2Is a 3X 3 size matrix, each column represents a coordinate axis unit vector of the coordinate system of the X-ray image captured by F2; v2The image coordinate system is a 3 multiplied by 1 size column vector, represents the original point position of the F2 image coordinate system, and takes the image center as the original point position of the F2 image coordinate system;
s124: establishing a rigid orthogonal space coordinate system (R) by using any lead point position on the positioner as an origin0，V0) Wherein R is0Is a 3 x 3 size matrix, each column representing a coordinate axis unit vector of the locator coordinate system; v0The coordinate system is a 3 multiplied by 1 row vector and represents the origin position of the coordinate system of the positioner;
s125: the position of each lead point on the locator can be represented by a column vector, and the positions of the 4 lead points forming a square on the locator relative to the locator coordinate system can be represented by a matrix P0Is represented by P0Is a 3 x 4 size matrix, each column representing 1 lead point coordinate;
relative to the positions of F1, F2, P01And P02Can be represented by a locator coordinate system (R)0，V0) To the F1 coordinate system (R)1，V1) And F2 coordinate system (R)2，V2) The coordinate system transformation between the two is expressed as follows:
P01＝R1 -1(R0P0+V0-V1)＝R01P0+V01；
P02＝R2 -1(R0P0+V0-V2)＝R02P0+V02；
wherein, P01Is a matrix representing the position of the lead points relative to F1, P02Is a matrix representing the position of the lead dots relative to F2;
s126: in the F1 image coordinate system, there is a projection matrix T based on the relative position of the F1 radioactive source and the F1 projection plane1，T1Is a 3 × 4 size matrix, P01Can be projected on an F1 plane to obtain F1 plane coordinates P01_im，P01_imIn a 2 x 4 size matrix, each column represents the coordinates of 1 lead point in the F1 plane as shown in the following equation:
extracting lead point coordinate P in F1 imageim1Calculating the projection coordinate P of the lead point at F1 by using a quasi-Newton method01_imAnd Pim1Coordinate system transformation (R) of locator-F1 when deviation between the two reaches minimum01，V01) As shown in the following formula:
s127: similarly, in the F2 image coordinate system, there is a projection matrix T based on the relative position of the F2 radioactive source and the F2 projection plane2，T2Is a 3 × 4 size matrix, P02Can be projected on an F2 plane to obtain F2 plane coordinates P02_im，P02_imIn a 2 x 4 size matrix, each column represents the coordinates of 1 lead point in the F2 plane as shown in the following equation:
extracting lead point coordinate P in F2 imageim2Calculating the projection coordinate P of the lead point at F2 by using a quasi-Newton method02_imAnd the lead point coordinate P in the real F2im2Coordinate system transformation (R) of locator-F2 when deviation between the two reaches minimum02，V02) As shown in the following formula:
F1-F2 coordinate system transformation (R)21，V21) Can be represented by the following formula:
R21＝R01R02 -1，V21＝V01-R01R02 -1V02；
s128: the spatial relationship between the two X-ray machine image planes can be used (R)21，V21) Showing based on the relative position of the source of each X-ray machine to the image plane, and the relative positional relationship (R) of the two X-ray machines21，V21) A virtual biplane X-ray system for automatic registration may be constructed in a computer system;
in the virtual biplane X-ray system, if the F1 image coordinate system is regarded as the global coordinate system for convenience of representation, the F1 image coordinate system (R) is regarded as the global coordinate system1，V1) And F2 image coordinate system (R)2，V2) Comprises the following steps:
R2＝R21，V2＝V21。
optionally, the step S2 further includes:
s21: taking an X-ray image shot by F1 and an X-ray image shot by F2 recorded at the same time as a group, and sequentially guiding each group of X-ray images into an F1 image plane and an F2 image plane; the X-ray image shot by F1 and the X-ray image shot by F2 are both 1024 × 1024 gray scale images, the gray scale value ranges from 0 to 255, and the gray scale values can be represented by 1024 × 1024 size matrixes;
s22: extracting F1 and F2 respectively shot in the current group by using Canny operatorImplant edges in X-ray images; canny operator parameter selection Tl＝0.01，Th0.15, and sigma is not less than 5;
the extraction result is a binary matrix with the same size of 1024 multiplied by 1024, wherein the edge is 1, and the non-edge is 0;
s23: taking the edge of the implant extracted by the Canny operator as an initial image I0Performing the circulation of the formula (1) to complete the binary image expansion; after k times of dilation, the initial edge image I is processed0And the generated dilated image IiSumming the images, normalizing to obtain a weighted expansion edge matrix B for subsequent similarity evaluation; the number of swelling times k is 5;
wherein:
obtaining a weighted expansion edge matrix B for an X-ray image taken at F11(ii) a Obtaining a weighted expansion edge matrix B for an X-ray image taken at F22。
Optionally, the step S3 further includes:
s31: introducing an implant model: the implant model includes an implant orthogonal coordinate system (R)obj，Vobj) Set of nodes PobjAnd a set of triangular patches Fobj(ii) a Wherein:
set of nodes PobjIs an N x 3 size matrix, each row represents 1 model node in the orthogonal coordinate system (R) of the implantobj，Vobj) The lower space coordinate, the number of rows N is the total number of nodes;
set of triangular patches FobjThe method is characterized in that the method is an M multiplied by 3 matrix, each line represents the serial number of nodes forming 1 triangular patch, and the line number M is the total number of patches;
adjusting an implant orthogonal coordinate system (R)obj，Vobj) The spatial positions of all implant model nodes in the virtual biplane X-ray system can be adjusted;
s32: based on an orthogonal coordinate system (R) of the implantobj，Vobj) And F1 space coordinate system (R)1，V1) Determining a set of implant model nodes PobjPosition P relative to the F1 coordinate systemobj1Projection matrix T determined based on the relative positions of the F1 radiation source and the F1 projection plane1The model node can be projected to the F1 plane to obtain the coordinate position of the model node in the plane
As shown in the following formula:
Pobj1＝R1 -1(RobjPobj+Vobj-V1)
s33: generating a 1024X 1024 zero matrix, spatially overlapping the X-ray image taken at F1, and collecting F according to whether the center is at the triangular patchobjAssigning 1 or 0 to each matrix element in the formed union set, wherein 1 represents that the center is in the triangular patch set, and 0 represents that the center is not in the triangular patch set, and generating a projection matrix Ap1And extracting an edge Ep1As shown in the following formula:
s34: similarly, based on the orthogonal coordinate system (R) of the implantobj，Vobj) And F2 space coordinate system (R)2，V2) A set of implant model nodes P may be determinedobjPosition P relative to the F2 coordinate systemobj2Projection matrix T determined based on the relative positions of the F2 radiation source and the F2 projection plane2The model node can be projected to the F2 plane to obtain the coordinate position of the model node in the plane
As shown in the following formula:
Pobj2＝R2 -1(RobjPobj+Vobj-V2)
s35: generating a 1024X 1024 zero matrix, spatially overlapping the X-ray image taken at F2, and collecting F according to whether the center is at the triangular patchobjAssigning 1 or 0 to each matrix element in the formed union set, wherein 1 represents that the center is in the triangular patch set, and 0 represents that the center is not in the triangular patch set, and generating a projection matrix Ap2And extracting an edge Ep2As shown in the following formula:
optionally, the step S4 further includes:
s41: in a constructed virtual biplane X-ray system, the operator manually adjusts the implant orthogonal coordinate system (R)obj，Vobj) To adjust the spatial position of the implant model nodes:
taking the F1 radioactive source as a camera position, the projection of the implant model on the F1 image plane is close to the implant area in the X-ray image shot by the F1; taking the F2 radioactive source as a camera position, the projection of the implant model on the F2 image plane is close to the implant area in the X-ray image shot by the F2; (ii) a Orthogonal coordinate system (R) of implant at this timeobj，Vobj) As an initial value for biplane auto-registration;
s42: computing a contour feature based matching function f (R)obj，Vobj) As shown in the following formula:
s43: when the degree of matching between the X-ray images captured by the F1 and the F2 and the projections of the models on the F1 and F2 image planes is maximum, the matching function F (R)obj，Vobj) Reaching a minimum value;
estimation of current model position (R) using simulated annealingobj，Vobj) Matching function f (R)obj，Vobj) Whether the minimum value is reached: when matching function f (R)obj，Vobj) If the minimum value is not reached, calculating iteration step length according to the annealing temperature so as to update the position of the model, and calculating a matching function value again according to S42; when matching function f (R)obj，Vobj) The position of the model at which the minimum value is reached (R)obj，Vobj) The real space position of the model in the virtual biplane X-ray system is obtained;
setting an initial temperature T020, the temperature reduction function selects an exponential-descent type, and the shutdown condition is set to Δ f < 1e-6(ii) a An objective function:
compared with the prior art, the invention has the following beneficial effects:
the invention popularizes and applies the 2D-3D automatic registration positioning method based on the single-plane X-ray to the biplane X-ray system, realizes the automatic registration positioning program of the implant based on the biplane X-ray tracking, and reduces the manpower and time cost of the manual registration of the biplane X-ray system.
Drawings
Other features, objects and advantages of the invention will become more apparent upon reading of the detailed description of non-limiting embodiments with reference to the following drawings:
FIG. 1 is a flow chart of a method for automatically registering and positioning an implant based on biplane X-ray tracking according to an embodiment of the present invention;
FIG. 2 is a schematic diagram of a virtual biplane X-ray system constructed in accordance with an embodiment of the present invention;
FIG. 3 is a schematic diagram of a virtual biplane X-ray system constructed in accordance with an embodiment of the present invention;
FIG. 4 is a schematic diagram of virtual biplane automatic registration according to an embodiment of the present invention;
FIG. 5 is a diagram illustrating a similarity evaluation method according to an embodiment of the present invention.
Detailed Description
The present invention will be described in detail with reference to specific examples. The following examples will assist those skilled in the art in further understanding the invention, but are not intended to limit the invention in any way. It should be noted that it would be obvious to those skilled in the art that various changes and modifications can be made without departing from the spirit of the invention. All falling within the scope of the present invention.
Referring to fig. 1 to 5, the present embodiment discloses an implant automatic registration positioning method based on biplane X-ray tracking, which includes the following steps:
s1: constructing a virtual biplane X-ray system, comprising:
s11: preparing two X-ray machines: the X-ray machine comprises a No. 1X-ray machine and a No. 2X-ray machine, and is used for shooting X-ray images of implants in a patient body;
s12: in a computer system, determining the relative positions of two X-ray machines;
s13: placing an implant model in a computer system in a public area projected by two X-ray machines;
s2: extracting implant edge features in the X-ray image of the implant, and performing weighted expansion;
s3: carrying out rapid projection and edge feature extraction on the implant model;
s4: and performing virtual biplane automatic registration based on the edge features obtained in the steps S2 and S3.
Wherein the S11 further comprises:
two X-ray machines were placed in the laboratory space: no. 1X-ray machine, No. 2X-ray machine make the patient who accomplishes the implant implantation carry out the function motion under the operator guide, and two X-ray machines shoot the implant and implant the joint simultaneously, obtain the X image respectively.
Wherein the step S12 further includes:
s121: a positioner (an acrylic plate is adopted in the embodiment) is arranged in a projection common area of two X-ray machines, the size of the acrylic plate is 100 multiplied by 100mm, 4 lead points are arranged in the middle of the positioner and distributed to form a 75 multiplied by 75mm square, the lead points have X-ray absorption property, and the lead points can be represented by low-gray circular points in an X-ray picture;
s122: the number 1X-ray machine is marked as F1, and a rigid orthogonal space coordinate system (R) is established by taking an X-ray image space shot by F1 as a reference1，V1) Wherein R is1Is a 3X 3 size matrix, each column represents a coordinate axis unit vector of the coordinate system of the X-ray image captured by F1; v1The image coordinate system is a 3 multiplied by 1 size column vector, represents the original point position of the F1 image coordinate system, and takes the image center as the original point position of the F1 image coordinate system; when the F1 image space is taken as the system coordinate system, it can be known that:
s123: the No. 2X-ray machine is marked as F2, and a rigid orthogonal space coordinate system (R) is established by taking an X-ray image space shot by F2 as a reference2，V2) Wherein R is2Is a 3X 3 size matrix, each column represents a coordinate axis unit vector of the coordinate system of the X-ray image captured by F2; v2The image coordinate system is a 3 multiplied by 1 size column vector, represents the original point position of the F2 image coordinate system, and takes the image center as the original point position of the F2 image coordinate system;
s124: establishing a rigid orthogonal space coordinate system (R) by using any lead point position on the positioner as an origin0，V0) Wherein R is0Is a 3 x 3 size matrix, each column representing a coordinate axis unit vector of the locator coordinate system; v0The coordinate system is a 3 multiplied by 1 row vector and represents the origin position of the coordinate system of the positioner;
s125: the position of each lead point on the locator can be represented by a column vector, and the positions of the 4 lead points forming a square on the locator relative to the locator coordinate system can be represented by a matrix P0Is represented by P0Is a 3 x 4 size matrix, each column representing 1 lead point coordinate;
relative to the positions of F1, F2, P01And P02Can be represented by a locator coordinate system (R)0，V0) To the F1 coordinate system (R)1，V1) And F2 coordinate system (R)2，V2) The coordinate system transformation between the two is expressed as follows:
P01＝R1 -1(R0P0+V0-V1)＝R01P0+V01；
P02＝R2 -1(R0P0+V0-V2)＝R02P0+V02；
wherein, P01Is a matrix representing the position of the lead points relative to F1, P02Is a matrix representing the position of the lead dots relative to F2;
s126: in the F1 image coordinate system, there is a projection matrix T based on the relative position of the F1 radioactive source and the F1 projection plane1，T1Is a 3 × 4 size matrix, P01Can be projected on an F1 plane to obtain F1 plane coordinates P01_im，P01_imIn a 2 x 4 size matrix, each column represents the coordinates of 1 lead point in the F1 plane as shown in the following equation:
the above formula is subjected to homogeneous coordinate transformation, the left side of the formula is a homogeneous coordinate, and w is the format writing method of the homogeneous equation here, represents a transformation relation and does not represent a specific matrix.
Extracting lead point coordinate P in F1 imageim1Calculating the projection coordinate P of the lead point at F1 by using a quasi-Newton method01_imAnd Pim1Coordinate system transformation (R) of locator-F1 when deviation between the two reaches minimum01，V01) As shown in the following formula:
s127: similarly, in the F2 image coordinate system, there is a projection matrix T based on the relative position of the F2 radioactive source and the F2 projection plane2，T2Is a 3 × 4 size matrix, P02Can be projected on an F2 plane to obtain F2 plane coordinates P02_im，P02_imIn a 2 x 4 size matrix, each column represents the coordinates of 1 lead point in the F2 plane as shown in the following equation:
extracting lead point coordinate P in F2 imageim2Calculating the projection coordinate P of the lead point at F2 by using a quasi-Newton method02_imAnd the lead point coordinate P in the real F2im2Coordinate system transformation (R) of locator-F2 when deviation between the two reaches minimum02，V02) As shown in the following formula:
F1-F2 coordinate system transformation (R)21，V21) Can be represented by the following formula:
R21＝R01R02 -1，V21＝V01-R01R02 -1V02；
s128: the spatial relationship between the two X-ray machine image planes can be used (R)21，V21) Showing based on the relative position of the source of each X-ray machine to the image plane, and the relative positional relationship (R) of the two X-ray machines21，V21) A virtual biplane X-ray system for automatic registration may be constructed in a computer system;
in the virtual biplane X-ray system, if the F1 image coordinate system is regarded as the global coordinate system for convenience of representation, the F1 image coordinate system (R) is regarded as the global coordinate system1，V1) And F2 image coordinate system (R)2，V2) Comprises the following steps:
R2＝R21，V2＝V21。
wherein the step S2 further includes:
s21: taking an X-ray image shot by F1 and an X-ray image shot by F2 recorded at the same time as a group, and sequentially guiding each group of X-ray images into an F1 image plane and an F2 image plane; the X-ray image shot by F1 and the X-ray image shot by F2 are both 1024 × 1024 gray scale images, the gray scale value ranges from 0 to 255, and the gray scale values can be represented by 1024 × 1024 size matrixes;
s22: extracting the edge of the implant in the X-ray images respectively shot by F1 and F2 in the current group by using a Canny operator; canny operator parameter selection Tl＝0.01，Th0.15, and sigma is not less than 5; in the present embodiment, σ is 3.
The parameter sigma controls the denoising effect and the blurring effect of the image, and sigma generally does not exceed 5 in order to keep the real implant edge information; the extraction result is a binary matrix with the same size of 1024 multiplied by 1024, wherein the edge is 1, and the non-edge is 0;
s23: to improve the robustness of the similarity assessment method, a weighted dilated F1, F2 edge image is used as input for the similarity assessment. Taking the edge of the implant extracted by the Canny operator as an initial image I0Performing the circulation of the formula (1) to complete the binary image expansion; after k times of dilation, the initial edge image I is processed0And the generated dilated image IiSumming the images, normalizing to obtain a weighted expansion edge matrix B for subsequent similarity evaluation; the number of swelling times k is 5; wherein:
obtaining a weighted expansion edge matrix B for an X-ray image taken at F11(ii) a Obtaining a weighted expansion edge matrix B for an X-ray image taken at F22。
Wherein the step S3 further includes:
s31: introducing an implant model: the implant model includes an implant orthogonal coordinate system (R)obj，Vobj) Set of nodes PobjAnd a set of triangular patches Fobj(ii) a Wherein:
set of nodes PobjIs an N x 3 size matrix, each row represents 1 model node in the orthogonal coordinate system (R) of the implantobj，Vobj) The lower space coordinate, the number of rows N is the total number of nodes;
set of triangular patches FobjThe method is characterized in that the method is an M multiplied by 3 matrix, each line represents the serial number of nodes forming 1 triangular patch, and the line number M is the total number of patches;
adjusting an implant orthogonal coordinate system (R)obj，Vobj) The spatial positions of all implant model nodes in the virtual biplane X-ray system can be adjusted;
s32: based on an orthogonal coordinate system (R) of the implantobj，Vobj) And F1 space coordinate system (R)1，V1) Determining a set of implant model nodes PobjPosition P relative to the F1 coordinate systemobj1Projection matrix T determined based on the relative positions of the F1 radiation source and the F1 projection plane1The model node can be projected to the F1 plane to obtain the coordinate position of the model node in the plane
As shown in the following formula:
Pobj1＝R1 -1(RobjPobj+Vobj-V1)
s33: for model nodes in a plane, a triangular patch set F is formedobjThe node connection relationship of (1) is not changed. Generating a 1024X 1024 zero matrix, spatially overlapping the X-ray image taken at F1, and collecting F according to whether the center is at the triangular patchobjAssigning 1 or 0 to each matrix element in the formed union set, wherein 1 represents that the center is in the triangular patch set, and 0 represents that the center is not in the triangular patch set, and generating a projection matrix Ap1And extracting an edge Ep1As shown in the following formula:
s34: similarly, based on the orthogonal coordinate system (R) of the implantobj，Vobj) And F2 space coordinate system (R)2，V2) A set of implant model nodes P may be determinedobjPosition P relative to the F2 coordinate systemobj2The projection matrix T2 determined based on the relative positions of the F2 radioactive source and the F2 projection plane can project the model nodes to the F2 plane to obtain the coordinate positions of the model nodes in the plane
As shown in the following formula:
Pobj2＝R2 -1(RobjPobj+Vobj-V2)
s35: for model nodes in a plane, a triangular patch set F is formedobjThe node connection relationship of (1) is not changed. Generating a 1024X 1024 zero matrix, spatially overlapping the X-ray image taken at F2, and collecting F according to whether the center is at the triangular patchobjAssigning 1 or 0 to each matrix element in the formed union, wherein 1 represents that the center is in the set of triangular patches and 0 represents that the center is not in the triangular patchesWithin the slice set, a projection matrix A is generatedp2And extracting an edge Ep2As shown in the following formula:
wherein the step S4 further includes:
s41: in a constructed virtual biplane X-ray system, the operator manually adjusts the implant orthogonal coordinate system (R)obj，Vobj) To adjust the spatial position of the implant model nodes:
taking the F1 radioactive source as a camera position, the projection of the implant model on the F1 image plane is close to the implant area in the X-ray image shot by the F1; taking the F2 radioactive source as a camera position, the projection of the implant model on the F2 image plane is close to the implant area in the X-ray image shot by the F2; the "close" criterion can be observed by the operator's eye with a deviation of no more than 100 pixels (10% of the image resolution). Orthogonal coordinate system (R) of implant at this timeobj，Vobj) As an initial value for biplane auto-registration.
S42: computing a contour feature based matching function f (R)obj，Vobj) As shown in the following formula:
s43: when the degree of matching between the X-ray images captured by the F1 and the F2 and the projections of the models on the F1 and F2 image planes is maximum, the matching function F (R)obj，Vobj) Reaching a minimum value;
estimation of current model position (R) using simulated annealingobj，Vobj) Matching function f (R)obj，Vobj) Whether the minimum value is reached: when matching function f (R)obj，Vobj) If the minimum value is not reached, calculating iteration step length according to the annealing temperature so as to update the position of the model, and calculating a matching function value again according to S42; when matching function f (R)obj，Vobj) The minimum value is reached (the stop condition is satisfied), and the model position (R) at this timeobj，Vobj) I.e. the real spatial position of the model in the virtual biplane X-ray system.
Setting an initial temperature T020, the temperature reduction function selects an exponential-descent type, and the shutdown condition is set to Δ f < 1e-6(ii) a An objective function:
the foregoing description of specific embodiments of the present invention has been presented. It is to be understood that the present invention is not limited to the specific embodiments described above, and that various changes or modifications may be made by one skilled in the art within the scope of the appended claims without departing from the spirit of the invention. The embodiments and features of the embodiments of the present application may be combined with each other arbitrarily without conflict.
Claims (6)
1. An implant automatic registration positioning method based on biplane X-ray tracking is characterized by comprising the following steps:
s1: constructing a virtual biplane X-ray system, comprising:
s11: preparing two X-ray machines: the X-ray machine comprises a No. 1X-ray machine and a No. 2X-ray machine, and is used for shooting X-ray images of implants in a patient body;
s12: in a computer system, determining the relative positions of two X-ray machines;
s13: placing an implant model in a computer system in a public area projected by two X-ray machines;
s2: extracting the edge of the implant in the X-ray image of the implant, and performing weighted expansion;
s3: performing fast projection and edge extraction on the implant model;
s4: and performing virtual biplane automatic registration based on the edge features obtained in the steps S2 and S3.
2. The method of claim 1, wherein the S11 further comprises:
two X-ray machines were placed in the laboratory space: no. 1X-ray machine, No. 2X-ray machine make the patient who accomplishes the implant implantation carry out the function motion under the operator guide, and two X-ray machines shoot the implant and implant the joint simultaneously, obtain the X image respectively.
3. The method of claim 1, wherein the step S12 further comprises:
s121: a positioner is arranged in a common projection area of two X-ray machines, 4 lead points are arranged in the middle of the positioner and distributed to form a square, the lead points have X-ray absorption and can be represented by low-gray circular points in an X-ray picture;
s122: the number 1X-ray machine is marked as F1, and a rigid orthogonal space coordinate system (R) is established by taking an X-ray image space shot by F1 as a reference1，V1) Wherein R is1Is a 3X 3 size matrix, each column represents a coordinate axis unit vector of the coordinate system of the X-ray image captured by F1; v1The image coordinate system is a 3 multiplied by 1 size column vector, represents the original point position of the F1 image coordinate system, and takes the image center as the original point position of the F1 image coordinate system; when the F1 image space is taken as the system coordinate system, it can be known that:
s123: the No. 2X-ray machine is marked as F2, and a rigid orthogonal space coordinate system (R) is established by taking an X-ray image space shot by F2 as a reference2，V2) Wherein R is2Is a 3X 3 size matrix, each column represents a coordinate axis unit vector of the coordinate system of the X-ray image captured by F2; v2The image coordinate system is a 3 multiplied by 1 size column vector, represents the original point position of the F2 image coordinate system, and takes the image center as the original point position of the F2 image coordinate system;
s124: establishing a rigid orthogonal space coordinate system (R) by using any lead point position on the positioner as an origin0，V0) Wherein R is0Is a 3 x 3 size matrix, each column representing a coordinate axis unit vector of the locator coordinate system; v0Is a column vector of size 3 x 1,representing the origin position of a coordinate system of the locator;
s125: the position of each lead point on the locator can be represented by a column vector, and the positions of the 4 lead points forming a square on the locator relative to the locator coordinate system can be represented by a matrix P0Is represented by P0Is a 3 x 4 size matrix, each column representing 1 lead point coordinate;
relative to the positions of F1, F2, P01And P02Can be represented by a locator coordinate system (R)0，V0) To the F1 coordinate system (R)1，V1) And F2 coordinate system (R)2，V2) The coordinate system transformation between the two is expressed as follows:
P01＝R1 -1(R0P0+V0-V1)＝R01P0+V01；
P02＝R2 -1(R0P0+V0-V2)＝R02P0+V02；
wherein, P01Is a matrix representing the position of the lead points relative to F1, P02Is a matrix representing the position of the lead dots relative to F2;
s126: in the F1 image coordinate system, there is a projection matrix T based on the relative position of the F1 radioactive source and the F1 projection plane1，T1Is a 3 × 4 size matrix, P01Can be projected on an F1 plane to obtain F1 plane coordinates P01_im，P01_imIn a 2 x 4 size matrix, each column represents the coordinates of 1 lead point in the F1 plane as shown in the following equation:
extracting lead point coordinate P in F1 imageim1Calculating the projection coordinate P of the lead point at F1 by using a quasi-Newton method01_imAnd Pim1Coordinate system transformation (R) of locator-F1 when deviation between the two reaches minimum01，V01) As shown in the following formula:
s127: similarly, in the F2 image coordinate system, there is a projection matrix T based on the relative position of the F2 radioactive source and the F2 projection plane2，T2Is a 3 × 4 size matrix, P02Can be projected on an F2 plane to obtain F2 plane coordinates P02_im，P02_imIn a 2 x 4 size matrix, each column represents the coordinates of 1 lead point in the F2 plane as shown in the following equation:
extracting lead point coordinate P in F2 imageim2Calculating the projection coordinate P of the lead point at F2 by using a quasi-Newton method02_imAnd the lead point coordinate P in the real F2im2Coordinate system transformation (R) of locator-F2 when deviation between the two reaches minimum02，V02) As shown in the following formula:
F1-F2 coordinate system transformation (R)21，V21) Can be represented by the following formula:
R21＝R01R02 -1，V21＝V01-R01R02 -1V02；
s128: the spatial relationship between the two X-ray machine image planes can be used (R)21，V21) Showing based on the relative position of the source of each X-ray machine to the image plane, and the relative positional relationship (R) of the two X-ray machines21，V21) A virtual biplane X-ray system for automatic registration may be constructed in a computer system;
in the virtual biplane X-ray system, the F1 image coordinate system is recorded as a global coordinate system for convenience of representationF1 image coordinate system (R)1，V1) And F2 image coordinate system (R)2，V2) Comprises the following steps:
R2＝R21，V2＝V21。
4. the method of claim 3, wherein the step S2 further comprises:
s21: taking an X-ray image shot by F1 and an X-ray image shot by F2 recorded at the same time as a group, and sequentially guiding each group of X-ray images into an F1 image plane and an F2 image plane; the X-ray image shot by F1 and the X-ray image shot by F2 are both 1024 × 1024 gray scale images, the gray scale value ranges from 0 to 255, and the gray scale values can be represented by 1024 × 1024 size matrixes;
s22: extracting the edge of the implant in the X-ray images respectively shot by F1 and F2 in the current group by using a Canny operator; canny operator parameter selection Tl＝0.01，Th0.15, and sigma is not less than 5;
the extraction result is a binary matrix with the same size of 1024 multiplied by 1024, wherein the edge is 1, and the non-edge is 0;
s23: taking the edge of the implant extracted by the Canny operator as an initial image I0Performing the circulation of the formula (1) to complete the binary image expansion; after k times of dilation, the initial edge image I is processed0And the generated dilated image IiSumming the images, normalizing to obtain a weighted expansion edge matrix B for subsequent similarity evaluation; the number of swelling times k is 5; wherein:
obtaining a weighted expansion edge matrix B for an X-ray image taken at F11(ii) a Obtaining a weighted expansion edge matrix B for an X-ray image taken at F22。
5. The method of claim 4, wherein the step S3 further comprises:
s31: introducing an implant model: the implant model includes an implant orthogonal coordinate system (R)obj，Vobj) Set of nodes PobjAnd a set of triangular patches Fobj(ii) a Wherein:
set of nodes PobjIs an N x 3 size matrix, each row represents 1 model node in the orthogonal coordinate system (R) of the implantobj，Vobj) The lower space coordinate, the number of rows N is the total number of nodes;
set of triangular patches FobjThe method is characterized in that the method is an M multiplied by 3 matrix, each line represents the serial number of nodes forming 1 triangular patch, and the line number M is the total number of patches;
adjusting an implant orthogonal coordinate system (R)obj，Vobj) The spatial positions of all implant model nodes in the virtual biplane X-ray system can be adjusted;
s32: based on an orthogonal coordinate system (R) of the implantobj，Vobj) And F1 space coordinate system (R)1，V1) Determining a set of implant model nodes PobjPosition P relative to the F1 coordinate systemobj1Projection matrix T determined based on the relative positions of the F1 radiation source and the F1 projection plane1The model node can be projected to the F1 plane to obtain the coordinate position of the model node in the plane
As shown in the following formula:
Pobj1＝R1 -1(RobjPobj+Vobj-V1)
s33: generating a 1024X 1024 zero matrix, spatially overlapping the X-ray image taken at F1, and collecting F according to whether the center is at the triangular patchobjAssigning 1 or 0 to each matrix element in the formed union set, wherein 1 represents that the center is in the triangular patch set, and 0 represents that the center is not in the triangular patch set, and generating a projection matrix Ap1And extracting an edge Ep1As shown in the following formula:
s34: similarly, based on the orthogonal coordinate system (R) of the implantobj，Vobj) And F2 space coordinate system (R)2，V2) A set of implant model nodes P may be determinedobjPosition P relative to the F2 coordinate systemobj2Projection matrix T determined based on the relative positions of the F2 radiation source and the F2 projection plane2The model node can be projected to the F2 plane to obtain the coordinate position of the model node in the plane
As shown in the following formula:
Pobj2＝R2 -1(RobjPobj+Vobj-V2)
s35: generating a 1024X 1024 zero matrix, spatially overlapping the X-ray image taken at F2, and collecting F according to whether the center is at the triangular patchobjAssigning 1 or 0 to each matrix element in the formed union set, wherein 1 represents that the center is in the triangular patch set, and 0 represents that the center is not in the triangular patch set, and generating a projection matrix Ap2And extracting an edge Ep2Of the formulaShown in the figure:
6. the method of claim 5, wherein the step S4 further comprises:
s41: in a constructed virtual biplane X-ray system, the operator manually adjusts the implant orthogonal coordinate system (R)obj，Vobj) To adjust the spatial position of the implant model nodes:
taking the F1 radioactive source as a camera position, the projection of the implant model on the F1 image plane is close to the implant area in the X-ray image shot by the F1; taking the F2 radioactive source as a camera position, the projection of the implant model on the F2 image plane is close to the implant area in the X-ray image shot by the F2; (ii) a Orthogonal coordinate system (R) of implant at this timeobj，Vobj) As an initial value for biplane auto-registration;
s42: computing a contour feature based matching function f (R)obj，Vobj) As shown in the following formula:
s43: when the degree of matching between the X-ray images captured by the F1 and the F2 and the projections of the models on the F1 and F2 image planes is maximum, the matching function F (R)obj，Vobj) Reaching a minimum value;
estimation of current model position (R) using simulated annealingobj，Vobj) Matching function f (R)obj，Vobj) Whether the minimum value is reached: when matching function f (R)obj，Vobj) If the minimum value is not reached, calculating iteration step length according to the annealing temperature so as to update the position of the model, and calculating a matching function value again according to S42; when matching function f (R)obj，Vobj) The position of the model at which the minimum value is reached (R)obj，Vobj) The real space position of the model in the virtual biplane X-ray system is obtained;
setting an initial temperature T020, the temperature reduction function selects an exponential-descent type, and the shutdown condition is set to Δ f < 1e-6(ii) a An objective function:
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010836044.3A CN111968164A (en) | 2020-08-19 | 2020-08-19 | Automatic implant registration and positioning method based on biplane X-ray tracking |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010836044.3A CN111968164A (en) | 2020-08-19 | 2020-08-19 | Automatic implant registration and positioning method based on biplane X-ray tracking |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN111968164A true CN111968164A (en) | 2020-11-20 |
Family
ID=73389398
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202010836044.3A Pending CN111968164A (en) | 2020-08-19 | 2020-08-19 | Automatic implant registration and positioning method based on biplane X-ray tracking |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN111968164A (en) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2023227511A1 (en) * | 2022-05-23 | 2023-11-30 | Koninklijke Philips N.V. | Simulating x-ray from low dose ct |
Citations (10)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US4630203A (en) * | 1983-12-27 | 1986-12-16 | Thomas Szirtes | Contour radiography: a system for determining 3-dimensional contours of an object from its 2-dimensional images |
| US4875165A (en) * | 1987-11-27 | 1989-10-17 | University Of Chicago | Method for determination of 3-D structure in biplane angiography |
| JPH05165956A (en) * | 1991-12-11 | 1993-07-02 | Fujitsu Ltd | Image processing method |
| US20030050527A1 (en) * | 2001-05-04 | 2003-03-13 | Peter Fox | Apparatus and methods for delivery of transcranial magnetic stimulation |
| CN105852971A (en) * | 2016-05-04 | 2016-08-17 | 苏州点合医疗科技有限公司 | Registration navigation method based on skeleton three-dimensional point cloud |
| CN105931190A (en) * | 2016-06-14 | 2016-09-07 | 西北工业大学 | High-angular-resolution light filed obtaining device and image generation method |
| CN106447704A (en) * | 2016-10-13 | 2017-02-22 | 西北工业大学 | A visible light-infrared image registration method based on salient region features and edge degree |
| CN107545585A (en) * | 2016-06-28 | 2018-01-05 | 西门子保健有限责任公司 | Method and apparatus for registering first image data set and the second image data set |
| CN110944594A (en) * | 2017-06-19 | 2020-03-31 | 穆罕默德·R·马赫福兹 | Hip surgical navigation using fluoroscopy and tracking sensors |
| CN111429532A (en) * | 2020-04-30 | 2020-07-17 | 南京大学 | Method for improving camera calibration accuracy by utilizing multi-plane calibration plate |
-
2020
- 2020-08-19 CN CN202010836044.3A patent/CN111968164A/en active Pending
Patent Citations (10)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US4630203A (en) * | 1983-12-27 | 1986-12-16 | Thomas Szirtes | Contour radiography: a system for determining 3-dimensional contours of an object from its 2-dimensional images |
| US4875165A (en) * | 1987-11-27 | 1989-10-17 | University Of Chicago | Method for determination of 3-D structure in biplane angiography |
| JPH05165956A (en) * | 1991-12-11 | 1993-07-02 | Fujitsu Ltd | Image processing method |
| US20030050527A1 (en) * | 2001-05-04 | 2003-03-13 | Peter Fox | Apparatus and methods for delivery of transcranial magnetic stimulation |
| CN105852971A (en) * | 2016-05-04 | 2016-08-17 | 苏州点合医疗科技有限公司 | Registration navigation method based on skeleton three-dimensional point cloud |
| CN105931190A (en) * | 2016-06-14 | 2016-09-07 | 西北工业大学 | High-angular-resolution light filed obtaining device and image generation method |
| CN107545585A (en) * | 2016-06-28 | 2018-01-05 | 西门子保健有限责任公司 | Method and apparatus for registering first image data set and the second image data set |
| CN106447704A (en) * | 2016-10-13 | 2017-02-22 | 西北工业大学 | A visible light-infrared image registration method based on salient region features and edge degree |
| CN110944594A (en) * | 2017-06-19 | 2020-03-31 | 穆罕默德·R·马赫福兹 | Hip surgical navigation using fluoroscopy and tracking sensors |
| CN111429532A (en) * | 2020-04-30 | 2020-07-17 | 南京大学 | Method for improving camera calibration accuracy by utilizing multi-plane calibration plate |
Non-Patent Citations (5)
| Title |
|---|
| 张俊华;吕梁;郭飞;吴俊;张榆锋;汪源源;施心陵;: "基于轮廓匹配的正交双平面X线图像重建胸段脊柱3D模型", 航天医学与医学工程, no. 04, 15 August 2012 (2012-08-15) * |
| 张斌 等: ""基于双平面DR 的骨折股骨三维姿态估计"", 《仪器仪表学报》, vol. 34, no. 1, 31 January 2013 (2013-01-31) * |
| 张烨明: ""机器人辅助髓内钉锁孔图像校正和定位算法研究"", 《中国优秀硕士学位论文全文数据库 (信息科技辑)》, 15 December 2006 (2006-12-15) * |
| 戚文潇;赵松;韩晓;赵金忠;: "双平面图像-模型配准技术在肩关节运动学研究中的应用", 国际骨科学杂志, no. 04, 25 July 2018 (2018-07-25) * |
| 方良毅 等: ""螺旋CT多平面重建技术对肺小结节术前定位的指导作用"", 临床放射学杂志, vol. 37, no. 4, 30 April 2018 (2018-04-30) * |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2023227511A1 (en) * | 2022-05-23 | 2023-11-30 | Koninklijke Philips N.V. | Simulating x-ray from low dose ct |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| US9681956B2 (en) | Acquiring and utilizing kinematic information for patient-adapted implants, tools and surgical procedures | |
| JP6294397B2 (en) | Customized orthopedic implants and related methods and deformable joint templates | |
| Tsai et al. | A volumetric model‐based 2D to 3D registration method for measuring kinematics of natural knees with single‐plane fluoroscopy | |
| Yamazaki et al. | Improvement of depth position in 2-D/3-D registration of knee implants using single-plane fluoroscopy | |
| Cappello et al. | Soft tissue artifact compensation in knee kinematics by double anatomical landmark calibration: performance of a novel method during selected motor tasks | |
| Baka et al. | 2D–3D shape reconstruction of the distal femur from stereo X-ray imaging using statistical shape models | |
| AU2009219989B2 (en) | Customised surgical apparatus | |
| Muhit et al. | Image-assisted non-invasive and dynamic biomechanical analysis of human joints | |
| Haque et al. | Hierarchical model-based tracking of cervical vertebrae from dynamic biplane radiographs | |
| CN115153830A (en) | System and method for assisting in assessing bone or soft tissue deformities using three-dimensional image reconstruction in orthopedic surgery | |
| CN111968164A (en) | Automatic implant registration and positioning method based on biplane X-ray tracking | |
| Blair et al. | Methodology for measurement of in vivo tibiotalar kinematics after total ankle replacement using dual fluoroscopy | |
| Ghanavati et al. | Multi-slice to volume registration of ultrasound data to a statistical atlas of human pelvis | |
| Roth et al. | An automated optimization pipeline for clinical-grade computer-assisted planning of high tibial osteotomies under consideration of weight-bearing | |
| WO2022229816A1 (en) | 3d reconstruction of anatomical images | |
| Saadat et al. | An efficient image registration method for 3D post-operative analysis of total knee arthroplasty | |
| US20230071033A1 (en) | Method for obtaining a ct-like representation and virtual x-ray images in arbitrary views from a two-dimensional x-ray image | |
| Saadat et al. | A fast and robust framework for 3D/2D model to multi-frame fluoroscopy registration | |
| US20240185509A1 (en) | 3d reconstruction of anatomical images | |
| US20230027518A1 (en) | Systems and methods for using photogrammetry to create patient-specific guides for orthopedic surgery | |
| Saadat | Multi-modal Image Registration | |
| Hampali | 3D Shape Reconstruction of Knee Bones from Low Radiation X-ray Images Using Deep Learning | |
| EP4246451A1 (en) | Method for modelling a joint | |
| Arn | An automated pipeline for three-dimensional preoperative planning of high tibial osteotomies under consideration of weight-bearing | |
| Hossain et al. | Repeat validation of a method to measure in vivo three dimensional hip kinematics using computed tomography and fluoroscopy |
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 | ||
| TA01 | Transfer of patent application right | ||
| TA01 | Transfer of patent application right |
Effective date of registration: 20211223 Address after: 201210 room 405-2, 4th floor, building 9, No. 1206, Zhangjiang Road, China (Shanghai) pilot Free Trade Zone, Pudong New Area, Shanghai Applicant after: SHANGHAI TAOYING MEDICAL TECHNOLOGY CO.,LTD. Address before: 200240 No. 800, Dongchuan Road, Shanghai, Minhang District Applicant before: SHANGHAI JIAO TONG University |