An automatic registration and positioning method for implants based on dual-plane X-ray tracking

Abstract

The present invention provides an automatic registration and positioning method for implants based on dual-plane X-ray tracking, including construction of a virtual dual-plane X-ray system, implant edge extraction and weighted expansion in X-ray images, rapid projection and edge extraction of implant geometric models, and search for optimal matching positions. According to a locator, the relative positions of the image planes of two X-ray machines are determined, thereby establishing a virtual dual-plane X-ray system in a computer system; the weighted expansion edges extracted from the X-ray images and the edges of the implant model projection are taken by two X-ray machines, and the matching degree of the implant model at any spatial position is calculated. The spatial position of the implant when the matching degree is the highest can be quickly obtained through a simulated annealing algorithm. This method can realize automatic tracking of the implant's spatial position in the body based on the motion images continuously obtained by the two X-ray machines, and use this as a reference for implant kinematic research and joint contact analysis, greatly reducing the manpower and time overhead of manual registration.

Description

Technical Field

[0001] The present invention relates to the fields of in-body kinematic tracking of metal implants, X-ray image post-processing, and 2D-3D image registration, and in particular to an automatic registration and positioning method for implants based on dual-plane X-ray tracking. Background Art

[0002] Total knee arthroplasty (TKA) is a method for treating end-stage knee osteoarthritis in clinical orthopedics. TKA mainly removes the defective bone tissue of the distal femur and proximal tibia, and reconstructs the joint using artificial implants made of metal, ceramic, high polyethylene and other materials, thereby eliminating the patient's knee pain. With the development of implant design principles, materials, configurations and other aspects, patients can restore their motor function to a certain extent after TKA, improve their quality of life and ability to live independently. Before surgery, the doctor will conduct a comprehensive assessment of factors such as the degree of progressive knee osteoarthritis, the function and integrity of the cruciate ligaments, and the patient's age, so as to select the most suitable surgical procedure and knee implant for the patient and decide on the implant configuration to be used. Clinical reports show that 15% to 30% of patients are dissatisfied with the postoperative effect of TKA[1]. Currently, the effective service life of implants after TKA is more than 90% of patients within 10-15 years, but the remaining patients who need implant revision bring a burden of 27 billion US dollars to the US health system[2]. The main factors that lead to the need for secondary revision surgery of knee implants include infection, implant loosening, gender, and age [1-3]. During surgery, the doctor determines the implant placement based on the joint state in a non-weight-bearing state. However, after surgery, the patient needs to perform a lot of rehabilitation activities such as sitting, standing up, walking, and going up and down stairs. The spatial movement and contact position of the knee implant are different from the static implant position during surgery, and the normal anatomical structure of the knee joint cannot be fully restored. Therefore, the implant may become loose and fail, requiring secondary surgery. Therefore, in vivo dynamic functional assessment of the knee joint is helpful to determine the optimal implant position of the knee implant, improve the implant configuration, and develop a personalized postoperative rehabilitation plan.

[0003] Dual Fluoroscopic Imaging System (DFIS) is a technology that uses two orthogonal X-ray machines to synchronously record joint motion images and accurately measure the spatial position and motion function of joints[4-6]. Two mobile X-ray machines are used to synchronously record the kinematic imaging data of the joints. Computer tomography (CT) and magnetic resonance imaging (MR) are used to obtain data to establish a personalized 3D model of the bone and knee implant. Based on the principle of point light source projection, the 3D bone model is matched with the 2D medical image in the virtual dual-plane system to obtain the spatial relative position of the joint for subsequent kinematic and other parameter analysis. The displacement and rotation angle accuracy of the dynamic motion of the knee joint measured by dual-plane dynamic X-ray tracking technology are 0.2 mm and 0.2° respectively[5]. However, in the current virtual dual-plane system, the operator needs to manually adjust the spatial position of the 3D model to achieve the 2D-3D registration process. The manual registration process is subjective, and the matching results of different operators may be different. In addition, it takes a lot of time. A skilled operator can quickly adjust the position of the 3D model to the vicinity of the precise position so that the projection contour and the 2D image contour roughly match (taking 10-20 seconds). However, in order to obtain the precise position of the 3D model, it takes 3-5 minutes to finely adjust the position of a single model. The knee joint is composed of the femur and tibia. It takes the operator more than 5 minutes to align the model of a set of 2D images, and it takes 300-400 minutes to fully align a set of motion data with a longer period (30-40 sets of 2D images). Therefore, it is very valuable to develop an automatic implant registration and positioning program based on dual-plane X-ray tracking to improve the repeatability and objectivity of the registration results and significantly reduce the manpower and time expenditure in the registration process.

[0004] At present, there are 2D-3D automatic registration and positioning methods based on single-plane X-rays, such as the automatic registration method of 2D single-plane X-rays and 3D metal implant models based on contour feature matching by Mahfouz et al. [7], the automatic registration method of 2D single-plane X-rays and 3D digital tomography data (CT) based on digitally reconstructed radiographs (DRRs) by Tsai et al. [8], and the real-time automatic registration method of 2D single-plane X-rays and 3D metal implant models based on convolutional neural networks by Miao et al. [9]. The basic principle of these automatic registration methods is to generate a virtual 2D plane projection image based on the relative positions of the 3D target, radiation source, and projection plane in a virtual single-plane projection system, calculate the similarity between the virtual 2D plane projection image and the actual transmission recorded 2D single-plane X-ray image, and use an optimization algorithm to adjust the spatial position of the 3D target so that the similarity between the projected 2D image and the actual recorded 2D single-plane X-ray image is the highest. There are three key technical issues in the technical implementation:

[0005] (1) Method for generating virtual 2D plane projection image.

[0006] ① In the computer system, 3D bone and metal implant models are composed of a series of "nodes" at specific positions and "interconnections between nodes" to form a set of 3D triangular facets. According to the basic theorem of perspective projection, the nodes and triangular facets can be directly projected onto the projection plane to generate a 2D binary mask to represent the projection result of the 3D model [7, 9].

[0007] ② CT is a 3D volume data set composed of a series of volume data units with the same spatial size. It can reconstruct the anatomical structure and relative position of bones and metal implants without weight bearing. The DRRs method

[10] can be used to generate 2D virtual X-ray fluoroscopic projection images based on the relationship between virtual radiation sources, CT data, and projection planes [8].

[0008] (2) Similarity evaluation method.

[0009] ① Similarity evaluation method based on grayscale distribution and contour features [7]. The virtual 2D plane projection image obtained based on the 3D model projection is a binary image (with projection is marked as 1, without projection is marked as 0). The actual recorded 2D single-plane X-ray image is an 8-bit grayscale image (grayscale value range 0-255). The normalized grayscale matching degree can be calculated using the following formula:

[0010]

[0011] Where G(x, y) is the actual recorded 2D single plane X-ray image gradient matrix, and H(x, y) is the virtual 2D plane projection image contour binary matrix. The comprehensive matching degree can be calculated by weighted summation of grayscale matching degree and contour matching degree.

[0012] The geometric contour of the virtual 2D plane projection image is extracted and represented by a binary matrix (contour is recorded as 1, non-contour is recorded as 0), and the gradient value of the actual recorded 2D single plane X-ray image is calculated. The gradient value of the contour is high and the gradient value of the non-contour is low. The normalized contour matching degree can be calculated using the following formula:

[0013]

[0014] Where J(x, y) is the grayscale matrix of the actual recorded 2D single-plane X-ray image, and J(x, y) is the binary matrix of the virtual 2D plane projection image.

[0015] ② Evaluation method based on gradient matching and weighted edge matching score (WEMS) [8]. The virtual 2D plane projection image projected based on the DRR method and the actual recorded 2D single plane X-ray image are both grayscale images, and their gradient values ​​can be calculated respectively. The gradient value matching degree is calculated using the following formula:

[0016]

[0017]

[0018]

[0019] Among them I f1 (x, y) is the grayscale matrix of the actual recorded 2D single-plane X-ray image, I DRR (x, y) is the grayscale matrix of the virtual 2D plane projection image generated by DRR.

[0020] The Canny edge detection operator is used to extract the edge of the actual recorded 2D single-plane X-ray image and generate a binary image (edge ​​is recorded as 1, non-edge is recorded as 0). The binary image expansion operation is used to perform weighted expansion on the binary edge image. Similarly, the Canny edge detection operator is used to extract the edge of the virtual 2D plane projection image and generate a binary image (edge ​​is recorded as 1, non-edge is recorded as 0). According to the length of the edge, the binary edge image is weighted and assigned a value. The longer the edge, the higher the weight. The weighted edge matching degree (WEMS) is calculated using the following formula:

[0021] E reg (x, y) = L DRR (x,y)B f1 (x, y).

[0022]

[0023] Among them B f1 (x, y) is the actual recorded 2D single plane X-ray image edge weight matrix after expansion, L DRR (x, y) is the edge weight matrix of the virtual 2D plane projection image after length assignment.

[0024] (3) Optimization calculation method.

[0025] ①Simulated Annealing (SA)[7];

[0026] ②Genetic Algorithms (GA)[8];

[0027] ③Convolutional Neural Network (CNN)[9].

[0028] The above prior art has the following defects:

[0029] (1) Traditional dual-plane dynamic X-ray imaging systems require manual registration, which is subjective and may produce different matching results for different operators. In addition, it takes a lot of time. A skilled operator can quickly adjust the position of the 3D model to near the precise position so that the projection contour and the 2D image contour roughly match (taking 10-20 seconds). However, in order to obtain the precise position of the 3D model, it takes 3-5 minutes to fine-tune the position of a single model. The knee joint is composed of the femur and tibia, and it takes the operator more than 5 minutes to align the model of a set of 2D images. It takes 300-400 minutes to fully align a set of motion data with a longer period (30-40 sets of 2D images).

[0030] (2) The 2D-3D automatic registration and positioning method based on single-plane X-ray has defects. In essence, single-plane X-ray is a single-view imaging system. When the position of the 3D target object is adjusted in the direction of the projection plane, the virtual 2D projection image will be significantly shifted, which significantly affects the matching degree between the virtual 2D image and the actual 2D image. However, when the position of the 3D target object is adjusted in the "out-of-plane" direction, the virtual 2D projection image will be enlarged or reduced in size. The enlargement or reduction ratio is affected by the size of the 3D target object's "out-of-plane" displacement and the distance between the radiation source and the projection plane. When the radiation source is close to the projection plane, a small displacement of the 3D target object in the "out-of-plane" direction can lead to a large size change of the virtual 2D projection image, but the space between the radiation source and the projection plane is small, and the perspective space is greatly limited; on the contrary, when the radiation source is far from the projection plane, the perspective space is guaranteed, but the small displacement of the 3D target object in the "out-of-plane" direction has almost no 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 during the functional movement of the subject, the radiation source of the X-ray machine is usually kept about 1 meter away from the projection plane. This makes the similarity evaluation function insensitive to the displacement of the 3D target in the "out-of-plane" direction during the 2D-3D registration process. Therefore, there is a large error in the "out-of-plane" direction when performing 2D-3D registration based on single-plane images. Mahfouz et al.'s automatic registration method of 2D single-plane X-ray and 3D metal implant model based on contour feature matching has an error of 1.4 mm in the "out-of-plane" direction [7]; Tsai et al.'s automatic registration method of 2D single-plane X-ray and 3D digital tomography (CT) data based on digitally reconstructed radiographic images has a maximum error of 7.4 mm in the "out-of-plane" direction when performing tibia bone registration [8]. This is because when extracting the edge of the 2D single-plane X-ray image of the bone, the bone edge of the non-target area will be extracted, which interferes with the result of the similarity evaluation function.

[0031] The references to the above content are as follows: The citation method is [number], where number is the corresponding numerical number in front of the following reference.

[0032] 1. Ilana N. Ackerman, Megan A. Bohensky, Ella Zomer, et al., The projected burst of primary total knee and hip replacement for osteoarthritis in Australia to the year 2030. Bmc Musculoskeletal Disorders, 2019.20(1).

[0033] 2.Douglas A.Dennis,Richard D.Komistek,and Mohamed R.Mahfouz,In VivoFluoroscopic Analysis Of Fixed-Bearing Total Knee Replacements.ClinicalOrthopaedics&Related Research,2003.410:p.114-130.

[0034] 3.Francesco Castagnini,Alessandra Sudanese,Barbara Bordini,et al.,Total Knee Replacement in Young Patients:Survival and Causes of Revision in ARegistry Population.Journal of Arthroplasty,2017:p.S0883540317304916.

[0035] 4.Tsung-Yuan Tsai,MingHanLincoln Liow,Guoan Li,et al.,Bi-CruciateRetaining Total Knee Arthroplasty Does Not Restore Native TibiofemoralArticular Contact Kinematics During Gait.Journal of Orthopaedic Research,2019.37(9).

[0036] 5.Jeffrey Bingham and Guoan Li,An Optimized Image Matching Method forDetermining In-Vivo TKA Kinematics with a Dual-Orthogonal FluoroscopicImaging System.Journal of Biomechanical Engineering,2006.128(4):p.588.

[0037] 6.Christophe Sauret,Hélène Pillet,Wafa Skalli,and Morgan Sangeux,Onthe use of knee functional calibration to determine the medio-lateral axis ofthe femur in gait analysis:Comparison with EOS biplanar radiographs asreference.Gait&Posture,2016.50.

[0038] 7.Mohamed R.Mahfouz,William A.Hoff,Richard D.Komistek,and DouglasA.Dennis,A robust method for registration of three-dimensional knee implantmodels to two-dimensional fluoroscopy images.IEEE Transactions on MedicalImaging,2003.22(12):p.p.1561-1574.

[0039] 8.Tsung-Yuan Tsai,Tung-Wu Lu,Chung-Ming Chen,et al.,A volumetricmodel-based 2D to 3D registration method for measuring kinematics of naturalknees with single-plane fluoroscopy.Medical Physics,2010.37(3).

[0040] 9.Shun Miao,Z.Jane Wang,and Liao Rui,A CNN Regression Approach forReal-Time 2D / 3D Registration.IEEE Transactions on Medical Imaging,2016.35(5):p.1352-1363.

[0041] 10.L.Lemieux, R.Jagoe, DRFish, et al., A patient-to-computed-tomographyimage registration method based on digitally reconstructed radiographs. Medical Physics, 1994.21(11): p.1749-1760. Summary of the invention

[0042] In view of the defects in the prior art, the purpose of the present invention is to provide an automatic registration and positioning method for implants based on dual-plane X-ray tracking. The technical solution of the present invention is as follows:

[0043] An automatic registration and positioning method for implants based on dual-plane X-ray tracking comprises the following steps:

[0044] S1: Build a virtual biplane X-ray system, including:

[0045] S11: preparing two X-ray machines: X-ray machine No. 1 and X-ray machine No. 2, and taking X-ray images of the implant in the patient's body;

[0046] S12: In the computer system, determining the relative positions of the two X-ray machines;

[0047] S13: Place the implant model in the common area projected by two X-ray machines in the computer system;

[0048] S2: extracting the implant edge in the X-ray image of the implant and performing weighted dilation;

[0049] S3: Rapid projection and edge extraction of implant models;

[0050] S4: Perform automatic virtual biplane registration based on the edge features obtained in steps S2 and S3.

[0051] Optionally, the S11 further includes:

[0052] Two X-ray machines are placed in the laboratory space: X-ray machine No. 1 and X-ray machine No. 2, so that patients who have completed implant placement can perform functional movements under the guidance of the operator, and at the same time the two X-ray machines take pictures of the implanted joints to obtain X-ray images.

[0053] Optionally, the step S12 further includes:

[0054] S121: Place a locator in the common area between the two X-ray machines. There are four lead dots in the middle of the locator, which are distributed in a square. The lead dots are X-ray absorbent and will be represented by low-grayscale dots in the X-ray image.

[0055] S122: X-ray machine No. 1 is denoted as F1, and a rigid body orthogonal space coordinate system (R1, V1) is established based on the X-ray image space taken by F1, where R1 is a 3×3 matrix, and each column represents the coordinate axis unit vector of the X-ray image coordinate system taken by F1; V1 is a 3×1 column vector, representing the origin position of the F1 image coordinate system, and the image center is the origin position of the F1 image coordinate system; when the F1 image space is used as the system coordinate system, it can be known that:

[0056]

[0057] S123: X-ray machine No. 2 is denoted as F2, and a rigid body orthogonal space coordinate system (R2, V2) is established based on the X-ray image space taken by F2, where R2 is a 3×3 matrix, and each column represents a unit vector of the coordinate axis of the X-ray image coordinate system taken by F2; V2 is a 3×1 column vector, representing the origin position of the F2 image coordinate system, and the image center is the origin position of the F2 image coordinate system;

[0058] S124: Taking any lead point position on the locator as the origin, establish a rigid body orthogonal space coordinate system (R0, V0), where R0 is a 3×3 matrix, each column of which represents the unit vector of the coordinate axis of the locator coordinate system; V0 is a 3×1 column vector, which represents the origin position of the locator coordinate system;

[0059] S125: The position of each lead point on the locator can be represented by a column vector, and the positions of the four lead points forming a square on the locator relative to the locator coordinate system can be represented by a matrix P0, where P0 is a 3×4 matrix, and each column represents the coordinate of one lead point;

[0060] Relative to the position of F1 and F2, P 01 and P 02 It can be expressed by the coordinate system transformation between the locator coordinate system (R0, V0) to the F1 coordinate system (R1, V1) and the F2 coordinate system (R2, V2), as shown in the following formula:

[0061] P 01 =R1 -1 (R0P0+V0-V1)=R 01 P0+V 01 ;

[0062] P 02 =R2 -1 (R0P0+V0-V2)=R 02 P0+V 02 ;

[0063] Among them, P 01 is the matrix representing the position of the lead point relative to F1, P02 is a matrix representing the position of the lead point relative to F2;

[0064] S126: In the F1 image coordinate system, based on the relative position of the F1 radiation source and the F1 projection plane, there is a projection matrix T1, T1 is a 3×4 matrix, P 01 It can be projected onto the F1 plane to obtain the F1 plane coordinate P 01_im , P 01_im It is a 2×4 matrix, each column represents the coordinates of a lead point on the F1 plane, as shown in the following formula:

[0065]

[0066] Extract the lead point coordinates P in the F1 image im1 , using the quasi-Newton method, calculate the projection coordinates P of the lead point on F1 01_im and P im1 The positioner-F1 coordinate system transformation (R 01 , V 01 ), as shown below:

[0067]

[0068] S127: Similarly, in the F2 image coordinate system, based on the relative position of the F2 radiation source and the F2 projection plane, there is a projection matrix T2, T2 is a 3×4 matrix, P 02 It can be projected onto the F2 plane to obtain the F2 plane coordinate P 02_im , P 02_im It is a 2×4 matrix, each column represents the coordinates of a lead point on the F2 plane, as shown in the following formula:

[0069]

[0070] Extract the lead point coordinates P in the F2 image im2 , using the quasi-Newton method, calculate the projection coordinates P of the lead point on F2 02_im and the lead point coordinate P in the real F2 im2 The positioner-F2 coordinate system transformation (R 02 , V 02 ), as shown below:

[0071]

[0072] F1-F2 coordinate transformation (R 21 , V 21 ) can be expressed by the following formula:

[0073] R 21 =R 01 R02 -1 , V 21 =V 01 -R 01 R 02 -1 V 02 ;

[0074] S128: The spatial relationship between the image planes of two X-ray machines can be expressed as (R 21 , V 21 ) indicates that the relative position of the radiation source of each X-ray machine and the image plane, as well as the relative position relationship of the two X-ray machines (R 21 , V 21 ), a virtual biplane X-ray system for automatic registration can be constructed in a computer system;

[0075] In the virtual dual-plane X-ray system, the F1 image coordinate system is recorded as the global coordinate system for the convenience of representation, and the F1 image coordinate system (R1, V1) and the F2 image coordinate system (R2, V2) are:

[0076]

[0077] R2=R 21 , V2=V 21 .

[0078] Optionally, step S2 further comprises:

[0079] S21: The X-ray images taken by F1 and the X-ray images taken by F2 recorded at the same time are grouped together, and each group of X-ray images is sequentially imported into the image planes of F1 and F2; wherein the X-ray images taken by F1 and the X-ray images taken by F2 are both 1024×1024 grayscale images, with grayscale values ​​ranging from 0 to 255, and can be represented by a 1024×1024 matrix;

[0080] S22: Use the Canny operator to extract the implant edges in the X-ray images taken by F1 and F2 in the current group; select T as the Canny operator parameter. l =0.01, T h =0.15, σ is not less than 5;

[0081] The extraction result is a binary matrix of the same size of 1024×1024, where the edge is 1 and the non-edge is 0;

[0082] S23: Using the implant edge extracted by the Canny operator as the initial image I0, perform a cycle of formula (1) to complete binary image expansion; after expansion k times, the initial edge image I0 and the generated expansion image I0 are combined. iThe images are summed and normalized to obtain the weighted dilation edge matrix B for subsequent similarity evaluation; here, the number of dilations k=5 is selected;

[0083] in:

[0084]

[0085]

[0086] For the X-ray image taken by F1, a weighted expansion edge matrix B1 is obtained; for the X-ray image taken by F2, a weighted expansion edge matrix B2 is obtained.

[0087] Optionally, step S3 further comprises:

[0088] S31: Import implant model: The implant model contains the implant orthogonal coordinate system (R obj , V obj ), node set P obj And the triangle patch set F obj ;in:

[0089] Node set P obj It is an N×3 matrix, each row represents a model node in the implant orthogonal coordinate system (R obj , V obj ) is the spatial coordinate under the node, and the number of rows N is the total number of nodes;

[0090] Triangle face set F obj It is an M×3 matrix, each row represents the serial number of the nodes that make up a triangular patch, and the number of rows M is the total number of patches;

[0091] Adjust the implant orthogonal coordinate system (R obj , V obj ), the spatial position of all implant model nodes in the virtual dual-plane X-ray system can be adjusted;

[0092] S32: Based on the implant orthogonal coordinate system (R obj , V obj ) and F1 space coordinate system (R1, V1) to determine the implant model node set P obj Position P relative to F1 coordinate system obj1 The projection matrix T1 determined based on the relative position of the F1 radiation source and the F1 projection plane can project the model nodes onto the F1 plane to obtain the coordinate position of the model nodes in the plane. As shown below:

[0093] P obj1 =R1 -1 (R obj P obj+V obj -V1)

[0094]

[0095] S33: Generate a zero matrix of size 1024×1024, which coincides with the X-ray image taken by F1 in spatial position, and calculates the value of the zero matrix according to whether the center is in the triangle face set F. obj Assign 1 or 0 to each matrix element in the union, where 1 indicates that the center is within the triangle set and 0 indicates that the center is not within the triangle set, generating the projection matrix A. p1 , and extract the edge E p1 , as shown below:

[0096]

[0097] S34: Similarly, based on the implant orthogonal coordinate system (R obj , V obj ) and F2 space coordinate system (R2, V2) can determine the implant model node set P obj Position P relative to F2 coordinate system obj2 The projection matrix T2 determined based on the relative position of the F2 radiation source and the F2 projection plane can project the model nodes onto the F2 plane to obtain the coordinate position of the model nodes in the plane. As shown below:

[0098] P obj2 =R2 -1 (R obj P obj +V obj -V2)

[0099]

[0100] S35: Generate a zero matrix of size 1024×1024, which coincides with the X-ray image taken by F2 in spatial position, and calculates the value of the zero matrix according to whether the center is in the triangle face set F. obj Assign 1 or 0 to each matrix element in the union, where 1 indicates that the center is within the triangle set and 0 indicates that the center is not within the triangle set, generating the projection matrix A. p2 , and extract the edge E p2 , as shown below:

[0101]

[0102] Optionally, the step S4 further comprises:

[0103] S41: In the 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:

[0104] With the F1 radiation source as the camera position, the projection of the implant model on the F1 image plane is close to the implant area in the X-ray image taken by F1; with the F2 radiation source as the camera position, the projection of the implant model on the F2 image plane is close to the implant area in the X-ray image taken by F2; at this time, the implant orthogonal coordinate system (R obj , V obj ) as the initial value for dual-plane automatic registration;

[0105] S42: Calculate the matching function f(R obj , V obj ), as shown below:

[0106]

[0107] S43: When the matching degree between the X-ray images taken by F1 and F2 and the projection of the model on the image planes of F1 and F2 reaches the maximum, the matching function f(R obj , V obj ) reaches the minimum value;

[0108] Use simulated annealing to evaluate the current model position (R obj , V obj )’s matching function f(R obj , V obj ) reaches the minimum value: When the matching function f(R obj , V obj ) does not reach the minimum value, the iteration step is calculated according to the annealing temperature at this time to update the model position, and the matching function value is recalculated according to S42; when the matching function f(R obj , V obj ) reaches its minimum value, and the model position (R obj , V obj ) is the real spatial position of the model in the virtual biplane X-ray system;

[0109] Set the initial temperature T0 = 20, select the exponential drop type for the cooling function, and set the shutdown condition to Δf < 1e -6 ; Objective function:

[0110]

[0111] Compared with the prior art, the present invention has the following beneficial effects:

[0112] The present invention promotes and applies the 2D-3D automatic registration and positioning method based on single-plane X-ray to a dual-plane X-ray system, realizes an implant automatic registration and positioning procedure based on dual-plane X-ray tracking, and reduces the manpower and time overhead of manual registration of the dual-plane X-ray system. BRIEF DESCRIPTION OF THE DRAWINGS

[0113] Other features, objects and advantages of the present invention will become more apparent from the detailed description of non-limiting embodiments made with reference to the following drawings:

[0114] Figure 1 It is a flow chart of an automatic registration and positioning method for implants based on dual-plane X-ray tracking according to a specific embodiment of the present invention;

[0115] Figure 2 is a schematic diagram of a virtual dual-plane X-ray system constructed according to a specific embodiment of the present invention;

[0116] Figure 3 is a schematic diagram of a virtual dual-plane X-ray system constructed according to a specific embodiment of the present invention;

[0117] Figure 4 is a schematic diagram of automatic virtual biplane registration according to a specific embodiment of the present invention;

[0118] Figure 5 It is a schematic diagram of a similarity evaluation method according to a specific embodiment of the present invention. DETAILED DESCRIPTION

[0119] The present invention is described in detail below in conjunction with specific embodiments. The following embodiments will help those skilled in the art to further understand the present invention, but are not intended to limit the present invention in any form. It should be noted that, for those of ordinary skill in the art, several changes and improvements can also be made without departing from the concept of the present invention. These all belong to the protection scope of the present invention.

[0120] See also Figures 1 to 5 This embodiment discloses an implant automatic registration and positioning method based on dual-plane X-ray tracking, comprising the following steps:

[0121] S1: Build a virtual biplane X-ray system, including:

[0122] S11: preparing two X-ray machines: X-ray machine No. 1 and X-ray machine No. 2, and taking X-ray images of the implant in the patient's body;

[0123] S12: In the computer system, determining the relative positions of the two X-ray machines;

[0124] S13: Place the implant model in the common area projected by two X-ray machines in the computer system;

[0125] S2: extracting implant edge features in the X-ray image of the implant and performing weighted dilation;

[0126] S3: fast projection and edge feature extraction of implant model;

[0127] S4: Perform automatic virtual biplane registration based on the edge features obtained in steps S2 and S3.

[0128] Wherein, the S11 further includes:

[0129] Two X-ray machines are placed in the laboratory space: X-ray machine No. 1 and X-ray machine No. 2, so that patients who have completed implant placement can perform functional movements under the guidance of the operator, and at the same time the two X-ray machines take pictures of the implanted joints to obtain X-ray images.

[0130] Wherein, the step S12 further comprises:

[0131] S121: A locator (an acrylic plate is used in this embodiment) is placed in the common projection area of ​​the two X-ray machines. The size of the acrylic plate is 100×100 mm. There are 4 lead dots in the middle of the locator, which are distributed to form a 75×75 mm square. The lead dots are X-ray absorbing and will be represented by low-grayscale dots in the X-ray image.

[0132] S122: X-ray machine No. 1 is denoted as F1, and a rigid body orthogonal space coordinate system (R1, V1) is established based on the X-ray image space taken by F1, where R1 is a 3×3 matrix, and each column represents the coordinate axis unit vector of the X-ray image coordinate system taken by F1; V1 is a 3×1 column vector, representing the origin position of the F1 image coordinate system, and the image center is the origin position of the F1 image coordinate system; when the F1 image space is used as the system coordinate system, it can be known that:

[0133]

[0134] S123: X-ray machine No. 2 is denoted as F2, and a rigid body orthogonal space coordinate system (R2, V2) is established based on the X-ray image space taken by F2, where R2 is a 3×3 matrix, and each column represents a unit vector of the coordinate axis of the X-ray image coordinate system taken by F2; V2 is a 3×1 column vector, representing the origin position of the F2 image coordinate system, and the image center is the origin position of the F2 image coordinate system;

[0135] S124: Taking any lead point position on the locator as the origin, establish a rigid body orthogonal space coordinate system (R0, V0), where R0 is a 3×3 matrix, each column of which represents the unit vector of the coordinate axis of the locator coordinate system; V0 is a 3×1 column vector, which represents the origin position of the locator coordinate system;

[0136] S125: The position of each lead point on the locator can be represented by a column vector, and the positions of the four lead points forming a square on the locator relative to the locator coordinate system can be represented by a matrix P0, where P0 is a 3×4 matrix, and each column represents the coordinate of one lead point;

[0137] Relative to the position of F1 and F2, P 01 and P 02 It can be expressed by the coordinate system transformation between the locator coordinate system (R0, V0) to the F1 coordinate system (R1, V1) and the F2 coordinate system (R2, V2), as shown in the following formula:

[0138] P 01 =R1 -1 (R0P0+V0-V1)=R 01 P0+V 01 ;

[0139] P 02 =R2 -1 (R0P0+V0-V2)=R 02 P0+V 02 ;

[0140] Among them, P 01 is the matrix representing the position of the lead point relative to F1, P 02 is a matrix representing the position of the lead point relative to F2;

[0141] S126: In the F1 image coordinate system, based on the relative position of the F1 radiation source and the F1 projection plane, there is a projection matrix T1, T1 is a 3×4 matrix, P 01 It can be projected onto the F1 plane to obtain the F1 plane coordinate P 01_im , P 01_im It is a 2×4 matrix, each column represents the coordinates of a lead point on the F1 plane, as shown in the following formula:

[0142]

[0143] The above formula performs a homogeneous coordinate transformation. The left side of the formula is a homogeneous coordinate. w here is the format of the homogeneous equation, which represents a transformation relationship, not a specific matrix.

[0144] Extract the lead point coordinates P in the F1 image im1 , using the quasi-Newton method, calculate the projection coordinates P of the lead point on F1 01_im and P im1 The positioner-F1 coordinate system transformation (R 01 , V 01 ), as shown below:

[0145]

[0146] S127: Similarly, in the F2 image coordinate system, based on the relative position of the F2 radiation source and the F2 projection plane, there is a projection matrix T2, T2 is a 3×4 matrix, P 02 It can be projected onto the F2 plane to obtain the F2 plane coordinate P 02_im , P 02_im It is a 2×4 matrix, each column represents the coordinates of a lead point on the F2 plane, as shown in the following formula:

[0147]

[0148] Extract the lead point coordinates P in the F2 image im2 , using the quasi-Newton method, calculate the projection coordinates P of the lead point on F2 02_im and the lead point coordinate P in the real F2 im2 The positioner-F2 coordinate system transformation (R 02 , V 02 ), as shown below:

[0149]

[0150] F1-F2 coordinate transformation (R 21 , V 21 ) can be expressed by the following formula:

[0151] R 21 =R 01 R 02 -1 , V 21 =V 01 -R 01 R 02 -1 V 02 ;

[0152] S128: The spatial relationship between the image planes of two X-ray machines can be expressed as (R 21 , V 21 ) indicates that the relative position of the radiation source of each X-ray machine and the image plane, as well as the relative position relationship of the two X-ray machines (R 21 , V 21 ), a virtual biplane X-ray system for automatic registration can be constructed in a computer system;

[0153] In the virtual dual-plane X-ray system, the F1 image coordinate system is recorded as the global coordinate system for the convenience of representation, and the F1 image coordinate system (R1, V1) and the F2 image coordinate system (R2, V2) are:

[0154]

[0155] R2=R 21 , V2=V 21 .

[0156] Wherein, the step S2 further comprises:

[0157] S21: The X-ray images taken by F1 and the X-ray images taken by F2 recorded at the same time are grouped together, and each group of X-ray images is sequentially imported into the image planes of F1 and F2; wherein the X-ray images taken by F1 and the X-ray images taken by F2 are both 1024×1024 grayscale images, with grayscale values ​​ranging from 0 to 255, and can be represented by a 1024×1024 matrix;

[0158] S22: Use the Canny operator to extract the implant edges in the X-ray images taken by F1 and F2 in the current group; select T as the Canny operator parameter. l =0.01, T h =0.15, σ is not less than 5; in this embodiment, σ=3.

[0159] The parameter σ controls the denoising and blurring effects of the image. In order to retain the real implant edge information, σ is generally not more than 5. The extraction result is a binary matrix of the same size of 1024×1024, where the edge is 1 and the non-edge is 0.

[0160] S23: In order to improve the robustness of the similarity evaluation method, the weighted expanded F1 and F2 edge images are used as the input of the similarity evaluation. The implant edge extracted by the Canny operator is used as the initial image I0, and the binary image expansion is completed by looping formula (1); after expansion k times, the initial edge image I0 and the generated expansion image I0 are combined. i The images are summed and normalized to obtain the weighted dilation edge matrix B for subsequent similarity evaluation; here, the number of dilations k = 5 is selected; where:

[0161]

[0162]

[0163] For the X-ray image taken by F1, a weighted expansion edge matrix B1 is obtained; for the X-ray image taken by F2, a weighted expansion edge matrix B2 is obtained.

[0164] Wherein, the step S3 further comprises:

[0165] S31: Import implant model: The implant model contains the implant orthogonal coordinate system (R obj , V obj ), node set P obj And the triangle patch set F obj ;in:

[0166] Node set P obj It is an N×3 matrix, each row represents a model node in the implant orthogonal coordinate system (R obj , V obj ) is the spatial coordinate under the node, and the number of rows N is the total number of nodes;

[0167] Triangle face set F obj It is an M×3 matrix, each row represents the serial number of the nodes that make up a triangular patch, and the number of rows M is the total number of patches;

[0168] Adjust the implant orthogonal coordinate system (R obj , V obj ), the spatial position of all implant model nodes in the virtual dual-plane X-ray system can be adjusted;

[0169] S32: Based on the implant orthogonal coordinate system (R obj , V obj ) and F1 space coordinate system (R1, V1) to determine the implant model node set P obj Position P relative to F1 coordinate system obj1 The projection matrix T1 determined based on the relative position of the F1 radiation source and the F1 projection plane can project the model nodes onto the F1 plane to obtain the coordinate position of the model nodes in the plane. As shown below:

[0170] P obj1 =R1 -1 (R obj P obj +V obj -V1)

[0171]

[0172] S33: For the model nodes in the plane, a triangular face set F is formed. obj The node connection relationship of has not changed. Generate a 1024×1024 size zero matrix, which coincides with the X-ray image taken by F1 in spatial position, and select the center according to whether it is in the triangle face set F obj Assign 1 or 0 to each matrix element in the union, where 1 indicates that the center is within the triangle set and 0 indicates that the center is not within the triangle set, generating the projection matrix A. p1 , and extract the edge E p1 , as shown below:

[0173]

[0174] S34: Similarly, based on the implant orthogonal coordinate system (R obj , Vobj ) and F2 space coordinate system (R2, V2) can determine the implant model node set P obj Position P relative to F2 coordinate system obj2 The projection matrix T2 determined based on the relative position of the F2 radiation source and the F2 projection plane can project the model nodes onto the F2 plane to obtain the coordinate position of the model nodes in the plane. As shown below:

[0175] P obj2 =R2 -1 (R obj P obj +V obj -V2)

[0176]

[0177] S35: For the model nodes in the plane, a triangular face set F is formed. obj The node connection relationship of has not changed. Generate a 1024×1024 size zero matrix, which coincides with the X-ray image taken by F2 in spatial position. obj Assign 1 or 0 to each matrix element in the union, where 1 indicates that the center is within the triangle set and 0 indicates that the center is not within the triangle set, generating the projection matrix A. p2 , and extract the edge E p2 , as shown below:

[0178]

[0179] Wherein, the step S4 further comprises:

[0180] S41: In the constructed virtual biplane X-ray system, the operator manually adjusts the implant orthogonal coordinate system (R obj , V obj ) to adjust the spatial position of the implant model nodes:

[0181] With the F1 radiation source as the camera position, the projection of the implant model on the F1 image plane is close to the implant area in the F1 X-ray image; with the F2 radiation source as the camera position, the projection of the implant model on the F2 image plane is close to the implant area in the F2 X-ray image; the above “close” standard can be based on the operator’s eye observation, and the deviation does not exceed 100 pixels (10% of the image resolution). At this time, the implant orthogonal coordinate system (R obj , V obj ) is used as the initial value for the automatic registration of the two planes.

[0182] S42: Calculate the matching function f(R obj, V obj ), as shown below:

[0183]

[0184] S43: When the matching degree between the X-ray images taken by F1 and F2 and the projection of the model on the image planes of F1 and F2 reaches the maximum, the matching function f(R obj , V obj ) reaches the minimum value;

[0185] Use simulated annealing to evaluate the current model position (R obj , V obj )’s matching function f(R obj , V obj ) reaches the minimum value: When the matching function f(R obj , V obj ) does not reach the minimum value, the iteration step is calculated according to the annealing temperature at this time to update the model position, and the matching function value is recalculated according to S42; when the matching function f(R obj , V obj ) reaches the minimum value (meets the stop condition), and the model position (R obj , V obj ) is the real spatial position of the model in the virtual dual-plane X-ray system.

[0186] Set the initial temperature T0 = 20, select the exponential drop type for the cooling function, and set the shutdown condition to Δf < 1e -6 ; Objective function:

[0187]

[0188] The above describes the specific embodiments of the present invention. It should be understood that the present invention is not limited to the above specific embodiments, and those skilled in the art can make various changes or modifications within the scope of the claims, which does not affect the essence of the present invention. In the absence of conflict, the embodiments of the present application and the features in the embodiments can be combined with each other arbitrarily.

Claims

1. An automatic registration and positioning method for implants based on dual-plane X-ray tracking, characterized in that: The steps include: S1: Build a virtual biplane X-ray system, including: S11: preparing two X-ray machines: X-ray machine No. 1 and X-ray machine No. 2, and taking X-ray images of the implant in the patient's body; S12: In the computer system, determining the relative positions of the two X-ray machines; S13: Place the implant model in the common area projected by two X-ray machines in the computer system; S2: extracting the implant edge in the X-ray image of the implant and performing weighted dilation; The step S2 further comprises: S21: record the X-ray machine No. 1 as F1, record the X-ray machine No. 2 as F2, take the X-ray images taken by F1 and the X-ray images taken by F2 recorded at the same time as a group, and import each group of X-ray images into the image planes of F1 and F2 in sequence; wherein the X-ray images taken by F1 and the X-ray images taken by F2 are both 1024×1024 grayscale images with a grayscale value range from 0 to 255, and can be represented by a 1024×1024 matrix; S22: Use the Canny operator to extract the implant edges in the X-ray images taken by F1 and F2 in the current group; select T as the Canny operator parameter. l =0.01,T h =0.15,σnot less than 5; The extraction result is a binary matrix of the same size of 1024×1024, where the edge is 1 and the non-edge is 0; S23: Using the implant edge extracted by the Canny operator as the initial image I0, perform a cycle of formula (1) to complete binary image expansion; after expansion k times, the initial edge image I0 and the generated expansion image I0 are combined. i The images are summed and normalized to obtain the weighted dilation edge matrix B for subsequent similarity evaluation; here, the number of dilations k = 5 is selected; where: For the X-ray image taken by F1, a weighted expansion edge matrix B1 is obtained; for the X-ray image taken by F2, a weighted expansion edge matrix B2 is obtained; S3: Rapid projection and edge extraction of implant models; S4: Perform automatic virtual biplane registration based on the edge features obtained in steps S2 and S3.

2. The method according to claim 1, characterized in that The S11 further comprises: Two X-ray machines are placed in the laboratory space: X-ray machine No. 1 and X-ray machine No. 2, so that patients who have completed implant placement can perform functional movements under the guidance of the operator, and at the same time the two X-ray machines take pictures of the implanted joints to obtain X-ray images.

3. The method according to claim 1, characterized in that The step S12 further comprises: S121: Place a locator in the common area between the two X-ray machines. There are four lead dots in the middle of the locator, which are distributed in a square. The lead dots are X-ray absorbent and will be represented by low-grayscale dots in the X-ray image. S122: A rigid body orthogonal space coordinate system (R1, V1) is established based on the X-ray image space taken by F1, where R1 is a 3×3 matrix, each column of which represents the coordinate axis unit vector of the X-ray image coordinate system taken by F1; V1 is a 3×1 column vector, which represents the origin position of the F1 image coordinate system, with the image center as the origin position of the F1 image coordinate system; when the F1 image space is used as the system coordinate system, it can be known that: S123: A rigid body orthogonal space coordinate system (R2, V2) is established based on the X-ray image space taken by F2, where R2 is a 3×3 matrix, each column of which represents the unit vector of the coordinate axis of the X-ray image coordinate system taken by F2; V2 is a 3×1 column vector, which represents the origin position of the F2 image coordinate system, with the image center as the origin position of the F2 image coordinate system; S124: Taking any lead point position on the locator as the origin, establish a rigid body orthogonal space coordinate system (R0, V0), where R0 is a 3×3 matrix, each column of which represents the unit vector of the coordinate axis of the locator coordinate system; V0 is a 3×1 column vector, which represents the origin position of the locator coordinate system; S125: The position of each lead point on the locator can be represented by a column vector, and the positions of the four lead points forming a square on the locator relative to the locator coordinate system can be represented by a matrix P0, where P0 is a 3×4 matrix, and each column represents the coordinate of one lead point; Relative to the position of F1 and F2, P 01 and P 02 It can be expressed by the coordinate system transformation between the locator coordinate system (R0, V0) to the F1 coordinate system (R1, V1) and the F2 coordinate system (R2, V2), as shown in the following formula: P 01 =R1 -1 (R0P0+V0-V1)=R 01 P0+V 01 ; <h2 style=";text-align:left;direction:ltr">P<h2 style=";text-align:left;direction:ltr"> 02 <h2 style=";text-align:left;direction:ltr"> =R2<h2 style=";text-align:left;direction:ltr"> -1 <h2 style=";text-align:left;direction:ltr"> (R0P0+V0-V2)=R<h2 style=";text-align:left;direction:ltr"> 02 <h2 style=";text-align:left;direction:ltr"> P0+V<h2 style=";text-align:left;direction:ltr"> 02 <h2 style=";text-align:left;direction:ltr"> ; Among them, P 01 is the matrix representing the position of the lead point relative to F1, P 02 is a matrix representing the position of the lead point relative to F2; S126: In the F1 image coordinate system, based on the relative position of the F1 radiation source and the F1 projection plane, there is a projection matrix T1, T1 is a 3×4 matrix, P 01 It can be projected onto the F1 plane to obtain the F1 plane coordinate P 01_im , P 01_im It is a 2×4 matrix, each column represents the coordinates of a lead point on the F1 plane, as shown in the following formula: Extract the lead point coordinates P in the F1 image im1 , using the quasi-Newton method, calculate the projection coordinates P of the lead point on F1 01_im and P im1 The positioner-F1 coordinate system transformation (R 01 ,V 01 ), as shown below: S127: Similarly, in the F2 image coordinate system, based on the relative position of the F2 radiation source and the F2 projection plane, there is a projection matrix T2, Y2 is a 3×4 size matrix, P 02 It can be projected onto the F2 plane to obtain the F2 plane coordinate P 02_im , P 02_im It is a 2×4 matrix, each column represents the coordinates of a lead point on the F2 plane, as shown in the following formula: Extract the lead point coordinates P in the F2 image im2 , using the quasi-Newton method, calculate the projection coordinates P of the lead point on F2 02_im and the lead point coordinate P in the real F2 im2 The positioner-F2 coordinate system transformation (R 02 ,V 02 ), as shown below: F1-F2 coordinate transformation (R 21 ,V 21 ) can be expressed by the following formula: R 21 =R 01 R 02 -1 ,V 21 =V 01 -R 01 R 02 -1 V 02 ; S128: The spatial relationship between the image planes of two X-ray machines can be expressed as (R 21 ,V 21 ) indicates that the relative position of the radiation source of each X-ray machine and the image plane, as well as the relative position relationship of the two X-ray machines (R 21 ,V 21 ), a virtual biplane X-ray system for automatic registration can be constructed in a computer system; In the virtual dual-plane X-ray system, the F1 image coordinate system is recorded as the global coordinate system for the convenience of representation, and the F1 image coordinate system (R1, V1) and the F2 image coordinate system (R2, V2) are: <h2 style=";text-align:left;direction:ltr">R2 = R<h2 style=";text-align:left;direction:ltr"> 21 <h2 style=";text-align:left;direction:ltr"> V2=V<h2 style=";text-align:left;direction:ltr"> 21 <h2 style=";text-align:left;direction:ltr"> 。 4. The method according to claim 3, characterized in that The step S3 further comprises: S31: Import implant model: The implant model contains the implant orthogonal coordinate system (R obj ,V obj ), node set P obj And the triangle patch set F obj ;in: Node set P obj It is an N×3 matrix, each row represents a model node in the implant orthogonal coordinate system (R obj ,V obj ) is the spatial coordinate under the node, and the number of rows N is the total number of nodes; Triangle face set F obj It is an M×3 matrix, each row represents the serial number of the nodes that make up a triangular patch, and the number of rows M is the total number of patches; Adjust the implant orthogonal coordinate system (R obj ,V obj ), the spatial position of all implant model nodes in the virtual dual-plane X-ray system can be adjusted; S32: Based on the implant orthogonal coordinate system (R obj ,V obj ) and F1 space coordinate system (R1, V1) to determine the implant model node set P obj Position P relative to F1 coordinate system obj1 The projection matrix T1 determined based on the relative position of the F1 radiation source and the F1 projection plane can project the model nodes onto the F1 plane to obtain the coordinate position of the model nodes in the plane. As shown below: P obj1 =R1 -1 (R obj P obj +V obj -V1) S33: Generate a zero matrix of size 1024×1024, which coincides with the X-ray image taken by F1 in spatial position, and calculates the value of the zero matrix according to whether the center is in the triangle face set F. obj Assign 1 or 0 to each matrix element in the union, where 1 indicates that the center is within the triangle set and 0 indicates that the center is not within the triangle set, generating the projection matrix A. p1 , and extract the edge E p1 , as shown below: E p1 =edge(A p1 ) S34: Similarly, based on the implant orthogonal coordinate system (R obj ,V obj ) and F2 space coordinate system (R2, V2) can determine the implant model node set P obj Position P relative to F2 coordinate system obj2 Based on the relative position of the F2 radiation source and the F2 projection plane, the projection matrix T2 can project the model nodes onto the F2 plane to obtain the coordinate position of the model nodes in the plane. As shown below: P obj2 =R2 -1 (R obj P obj +V obj -V2) S35: Generate a zero matrix of size 1024×1024, which coincides with the X-ray image taken by F2 in spatial position, and calculates the value of the zero matrix according to whether the center is in the triangle face set F. obj Assign 1 or 0 to each matrix element in the union, where 1 indicates that the center is within the triangle set and 0 indicates that the center is not within the triangle set, generating the projection matrix A. p2 , and extract the edge E p2 , as shown below: E p2 =edge(A p2 )。 5. The method according to claim 4, characterized in that The step S4 further comprises: S41: In the constructed virtual biplane X-ray system, the operator manually adjusts the implant orthogonal coordinate system (R obj ,V obj ) to adjust the spatial position of the implant model nodes: With the F1 radiation source as the camera position, the projection of the implant model on the F1 image plane is close to the implant area in the X-ray image taken by F1; with the F2 radiation source as the camera position, the projection of the implant model on the F2 image plane is close to the implant area in the X-ray image taken by F2; at this time, the implant orthogonal coordinate system (R obj ,V obj ) as the initial value for dual-plane automatic registration; S42: Calculate the matching function f(R obj ,V obj ), as shown below: S43: When the matching degree between the X-ray images taken by F1 and F2 and the projection of the model on the image planes of F1 and F2 reaches the maximum, the matching function f(R obj ,V obj ) reaches a minimum value; Use simulated annealing to evaluate the current model position (R obj ,V obj )’s matching function f(R obj ,V obj ) reaches the minimum value: When the matching function f(R obj ,V obj ) does not reach the minimum value, the iteration step is calculated according to the annealing temperature at this time to update the model position, and the matching function value is recalculated according to S42; when the matching function f(R obj ,V obj ) reaches its minimum value, and the model position (R obj ,V obj ) is the real spatial position of the model in the virtual biplane X-ray system; Set the initial temperature T0 = 20, select the exponential drop type for the cooling function, and set the shutdown condition to △f < 1e -6 ; Objective function:

Images