CN111568456A - Knee joint posture measuring method based on feature point three-dimensional reconstruction - Google Patents

Knee joint posture measuring method based on feature point three-dimensional reconstruction Download PDF

Info

Publication number
CN111568456A
CN111568456A CN202010331744.7A CN202010331744A CN111568456A CN 111568456 A CN111568456 A CN 111568456A CN 202010331744 A CN202010331744 A CN 202010331744A CN 111568456 A CN111568456 A CN 111568456A
Authority
CN
China
Prior art keywords
point
projection
knee joint
coordinate system
plane
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202010331744.7A
Other languages
Chinese (zh)
Other versions
CN111568456B (en
Inventor
师为礼
蒋振刚
李岩芳
许毅
苗语
赵家石
张科
何飞
秦俊
张丽媛
曲峰
何巍
杨晨
杨华民
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Changchun University of Science and Technology
Original Assignee
Changchun University of Science and Technology
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Changchun University of Science and Technology filed Critical Changchun University of Science and Technology
Priority to CN202010331744.7A priority Critical patent/CN111568456B/en
Publication of CN111568456A publication Critical patent/CN111568456A/en
Application granted granted Critical
Publication of CN111568456B publication Critical patent/CN111568456B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
    • A61B6/02—Arrangements for diagnosis sequentially in different planes; Stereoscopic radiation diagnosis
    • A61B6/03—Computed tomography [CT]
    • A61B6/032—Transmission computed tomography [CT]
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
    • A61B6/50—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment specially adapted for specific body parts; specially adapted for specific clinical applications
    • A61B6/505—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment specially adapted for specific body parts; specially adapted for specific clinical applications for diagnosis of bone
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects

Landscapes

  • Health & Medical Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Medical Informatics (AREA)
  • Physics & Mathematics (AREA)
  • Heart & Thoracic Surgery (AREA)
  • General Health & Medical Sciences (AREA)
  • High Energy & Nuclear Physics (AREA)
  • Theoretical Computer Science (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Optics & Photonics (AREA)
  • Pathology (AREA)
  • Radiology & Medical Imaging (AREA)
  • Biomedical Technology (AREA)
  • Veterinary Medicine (AREA)
  • Molecular Biology (AREA)
  • Surgery (AREA)
  • Animal Behavior & Ethology (AREA)
  • Biophysics (AREA)
  • Public Health (AREA)
  • Pulmonology (AREA)
  • Orthopedic Medicine & Surgery (AREA)
  • Dentistry (AREA)
  • Oral & Maxillofacial Surgery (AREA)
  • Computer Graphics (AREA)
  • Geometry (AREA)
  • Software Systems (AREA)
  • General Physics & Mathematics (AREA)
  • Image Analysis (AREA)
  • Image Processing (AREA)

Abstract

The invention relates to a knee joint posture measuring method based on feature point three-dimensional reconstruction, which comprises the following steps: establishing a three-dimensional model of the knee joint, and obtaining a spatial position parameter of an X-ray imaging system; extracting position information of projection points corresponding to the six characteristic points from the left X-ray image and the right X-ray image respectively; estimating the position coordinates of the feature points in the three-dimensional space; computing quaternion q from distal femur to projection pose1And quaternion q of proximal tibial transition to projection attitude2(ii) a From quaternion q1And quaternion q2And calculating to obtain a knee joint posture angle, and determining the posture of the knee joint. The invention avoids the influence of iterative registration process and soft tissue artifact by three-dimensionally reconstructing six characteristic points in the knee joint femur and tibia, has simple and convenient process, reduces calculated amount, and improvesHigh calculation efficiency and high accuracy.

Description

Knee joint posture measuring method based on feature point three-dimensional reconstruction
Technical Field
The invention relates to the technical field of medical auxiliary equipment, in particular to a knee joint posture measuring method based on feature point three-dimensional reconstruction.
Background
The knee joint is the most complex joint in the human body and the joint with the greatest degree of freedom, participates in most movements in daily life, and maintains the dynamic balance of the body. Because the knee joint structure is complicated, the amount of exercise is big, so the knee joint has very high sick and wounded rate in reality. The method is beneficial to more comprehensively understanding the relationship between the structure and the motion function of the knee joint for the basic research of the knee joint three-dimensional kinematics, and has important significance for researching the motion damage mechanism, clinical diagnosis, guiding rehabilitation training and the like. In the rehabilitation treatment stage, the kinematic parameters of the knee joint or the artificial joint and the lower limb movement posture can be detected non-invasively, and the method has very important significance for better monitoring and evaluating the rehabilitation condition of a patient by a doctor.
At present, the method for measuring the motion posture of the knee joint by using an optical tracking measuring device is a common method, but the method has low accuracy due to the influence of soft tissue artifacts. In addition, an in-vivo measurement technology utilizing 2D-3D registration is also a common method for measuring the motion posture of the knee joint, and although the method has higher accuracy than an optical tracking system, the process is complex and the calculation efficiency is relatively low.
Disclosure of Invention
The invention aims to provide a knee joint posture measuring method based on feature point three-dimensional reconstruction, aiming at the problems in the prior art, and the method has higher accuracy and calculation efficiency.
The invention adopts the following technical scheme:
a knee joint posture measuring method based on feature point three-dimensional reconstruction is characterized by comprising the following steps:
step one, carrying out CT scanning on the knee joint, and establishing a three-dimensional model of the knee joint; calibrating an X-ray imaging system with a left X-ray source, a right X-ray source and a corresponding imaging plane to obtain a spatial position parameter of the X-ray imaging system, and establishing a virtual imaging space in three-dimensional modeling software by using the obtained spatial position parameter;
selecting six characteristic points in the three-dimensional model, and extracting position information of projection points corresponding to the six characteristic points from a left X-ray image corresponding to the left X-ray source and a right X-ray image corresponding to the right X-ray source respectively;
thirdly, forming projection rays by using the projection points in the left X-ray source and the left X-ray image and the projection points in the right X-ray source and the right X-ray image, and estimating the position coordinates of the feature points in a three-dimensional space;
step four, forming three characteristic points on the distal femur in the three-dimensional model into a first plane, forming the three characteristic points into a second plane according to the estimated position coordinates of the three characteristic points on the distal femur, and aligning the three characteristic points with the second planeThe distal femur of the three-dimensional model is subjected to spatial transformation, so that the first plane is parallel to the second plane, and a quaternion q for transforming the distal femur into a projection posture is obtained through calculation1；
Forming three characteristic points on the proximal tibia end in the three-dimensional model into a third plane, forming the three characteristic points into a fourth plane according to estimated position coordinates of the three characteristic points on the proximal tibia end, performing spatial transformation on the proximal tibia end of the three-dimensional model to enable the third plane and the fourth plane to be parallel, and calculating to obtain a quaternion q converted from the proximal tibia end to a projection posture2；
Step five, according to the quaternion q1And said quaternion q2And calculating to obtain a knee joint posture angle, and determining the posture of the knee joint.
Compared with the prior art, the invention has the beneficial effects that: the knee joint posture measuring method based on the three-dimensional reconstruction of the characteristic points carries out the three-dimensional reconstruction of the six characteristic points in the femur and the tibia of the knee joint, avoids the influence of an iterative registration process and soft tissue artifacts, has simple and convenient process, reduces the calculated amount, improves the calculation efficiency and has higher accuracy.
Drawings
FIG. 1 is a flow chart of a knee joint posture measurement method based on feature point three-dimensional reconstruction according to the present invention;
FIG. 2 is a schematic spatial position diagram of an X-ray imaging system;
FIG. 3 is a schematic representation of feature points selected on a three-dimensional model of a knee joint.
Detailed Description
The technical solution of the present invention will be described in detail with reference to the accompanying drawings and preferred embodiments.
In one embodiment, as shown in fig. 1, the invention discloses a knee joint posture measurement method based on feature point three-dimensional reconstruction, which includes the following steps:
step one (S100), carrying out CT scanning on a knee joint (the distal end of a femur and the proximal end of a tibia) to establish a three-dimensional model of the knee joint; the X-ray imaging system comprises a left X-ray source, a right X-ray source and imaging planes corresponding to the left X-ray source and the right X-ray source, the X-ray imaging system is calibrated, namely the left X-ray source and the corresponding imaging planes as well as the right X-ray source and the corresponding imaging planes are calibrated, space position parameters of the X-ray imaging system are obtained, and a virtual imaging space is established in three-dimensional modeling software by utilizing the obtained space position parameters;
step two (S200), selecting six characteristic points (three characteristic points are respectively selected in the femur and the tibia) from the three-dimensional model of the knee joint, and extracting the position information of projection points corresponding to the six characteristic points from a left X-ray image corresponding to the left X-ray source and a right X-ray image corresponding to the right X-ray source respectively;
step three (S300), forming projection rays by using projection points in the left X-ray source and the left X-ray image and projection points in the right X-ray source and the right X-ray image, and estimating position coordinates of the characteristic points in a three-dimensional space;
step four (S400), three feature points on the distal femur in the three-dimensional model form a first plane, the three feature points form a second plane according to the estimated position coordinates of the three feature points on the distal femur, the distal femur of the three-dimensional model is subjected to space transformation to enable the first plane and the second plane to be parallel, and a quaternion q for transforming the distal femur into a projection posture is calculated1；
Forming three characteristic points on the proximal tibia end in the three-dimensional model into a third plane, forming the three characteristic points into a fourth plane according to estimated position coordinates of the three characteristic points on the proximal tibia end, carrying out spatial transformation on the proximal tibia end of the three-dimensional model to enable the third plane and the fourth plane to be parallel, and calculating to obtain a quaternion q converted from the proximal tibia end to a projection posture2；
Step five (S500), according to quaternion q1And quaternion q2And calculating to obtain a knee joint posture angle, and determining the posture of the knee joint.
As a specific implementation manner, in the step one (S100), the process of obtaining the spatial position parameter of the X-ray imaging system includes the following steps:
referring to FIG. 2, FIG. 2 shows an X-ray imaging system (i.e., two X-ray sources)Corresponding imaging plane), wherein point A represents a certain characteristic point on a three-dimensional model of the knee joint, points a1 and a2 represent projection points of the characteristic point A in the right imaging plane and the left imaging plane respectively, the left X-ray source and the right X-ray source are placed at an included angle of 30 degrees, the X-ray imaging system is subjected to stereo calibration by using a Zhang-friend calibration method, and an internal reference matrix K of the left X-ray source and the corresponding imaging plane is obtained1External reference matrix M1A camera matrix P1And an internal reference matrix K of the right X-ray source and the corresponding imaging plane2External reference matrix M2A camera matrix P2The positional relationship of the left and right X-ray sources may be represented by a matrix E as:
Figure BDA0002465188280000041
wherein x islIs the coordinate of a certain point X in space under the camera coordinate system of the left X-ray source, XrThe coordinate of a certain point X in the space under the camera coordinate system of the right X-ray source, R is a rotation matrix of the camera coordinate system of the right X-ray source and the camera coordinate system of the left X-ray source, and t is a translation vector;
after obtaining the spatial position parameters (internal reference matrix, external reference matrix and camera matrix) of the two X-ray sources and the imaging plane thereof, a basic matrix F of the X-ray imaging system can be obtained, wherein the basic matrix F is:
F＝[K2t]xK2R K1 -1(2)
wherein [ K ]2t]xIs an antisymmetric matrix of three-dimensional column vectors.
As a specific embodiment, the second step (S200) of selecting six feature points in the three-dimensional model of the knee joint and extracting the position information of the projection points corresponding to the six feature points in the left X-ray image and the right X-ray image respectively includes the following steps:
the selection of the knee joint feature points requires that three non-collinear mark points are respectively selected at the distal end of the femur and the proximal end of the tibia, the mark points should have obvious features, and the positions of projection points projected in an X-ray image are convenient to confirm, so that the medial epicondyle of the femur, the lateral epicondyle of the femur and the intercondylar notch of the femur on the distal end of the femur in the three-dimensional model are selected as three feature points of the femur, and the lateral condyle, the medial condyle and the intercondylar eminence of the tibia on the proximal end of the tibia in the three-dimensional model are selected as three;
manually marking the projection point positions of the characteristic points in the left X-ray image and the right X-ray image respectively;
and performing binarization processing on the left X-ray image and the right X-ray image after the projection point positions are marked, and respectively extracting the position information of the projection points.
As a specific implementation manner, after the position information of the projection points corresponding to the six feature points is extracted in step two (S200), the method further includes a process of calibrating the extracted position information of the projection points to obtain the position coordinates of the optimal projection point, where the process includes the following steps:
step two, assuming that the projection point extracted from the left X-ray image is X1＝(m，n，1)TThe projection point extracted in the right X-ray image is X2＝(m’，n’，1)TProjection point x1＝(m，n，1)TThe corresponding optimal projection point is x'1Projection point x2＝(m’，n’，1)TThe corresponding optimal projection point is x'2Projection point x1＝(m，n，1)TProjection point x2＝(m’，n’，1)TOptimal projected point x'1And optimal projected point x'2The distance minimization formula is satisfied:
C(x1，x2)＝d(x1、x’1)2+d(x’2，x’2)2(3)
step two, defining projection point x1＝(m，n，1)TAnd x2＝(m’，n’，1)TAnd translating it to the origin of coordinates:
Figure BDA0002465188280000051
Figure BDA0002465188280000052
step two and step three, respectively calculating a pole e in an imaging plane corresponding to the left X-ray source1＝(ex，ey，ez)TPole e in the imaging plane corresponding to the right X-ray source2＝(e’x，e’y，e’z)TPole e1And pole e2Respectively satisfy Fe10 and FTe20, pole e1And pole e2Respectively normalized to e2 x+e2 y＝1，e’x 2+e’y 21, and constructing a rotation matrix R1、R2So that it satisfies R1e1＝(1，0，ez)T、R2e1＝(1，0，e’z)TConstructed rotation matrix R1、R2Respectively as follows:
Figure BDA0002465188280000061
step two and four, the transformed basic matrix is F' ═ R2T2 -TF T1R1 TLet f be ez、f’＝e’z、a＝F’22、b＝F’23、c＝F’32、d＝F’33The transformed basis matrix F' may be represented in the form:
Figure BDA0002465188280000062
parameterizing epipolar line bundles in the left X-ray image by using a parameter u, wherein the epipolar line corresponding to the left X-ray image is set as l (u), and the epipolar line l (u) can be represented as a point (0, u, 1) passing through the left X-ray imageTWith transformed poles R in the left X-ray image1e1The cross product of (a), namely:
l(u)＝(0，u，1)T×(1，0，ez)T
＝(0，u，1)T×(1，0，f)T
＝(uf，1，-u)T(8)
let the epipolar line corresponding to the right X-ray image be l ' (u), and the epipolar line l ' (u) can be calculated by using the transformed fundamental matrix F ':
l’(u)＝F’(0，u，1)T＝(-f’(cu+d)，au+b，cu+d)T(9)
step two, because the projection point has already translated to the origin of coordinates in left X-ray image, right X-ray image, so projection point X1＝(m，n，1)TAnd projection point x2＝(m’，n’，1)TThe distance from the corresponding epipolar line in the left and right X-ray images, respectively, can be expressed as:
Figure BDA0002465188280000071
Figure BDA0002465188280000072
step two, because the optimal projection point is on the epipolar line, the formula (3) is expressed as a distance function parameterized by the parameter u as:
Figure BDA0002465188280000073
the derivation of equation (12) yields:
Figure BDA0002465188280000074
the denominators are combined and the numerator is made equal to 0 to obtain
g(u)＝u((au+b)2+f'2(cu+d)2)2-(1+f2u2)2(ad-bc)(au+b)(cu+d)＝0 (14)
Step two eight, solution type (14)Obtaining 6 roots, and comparing values of s (u) → ∞ time s (u) obtained by substituting the 6 roots into formula (12), thereby obtaining an optimal solution u (u) in which s (u) is minimizedmin；
Step two and nine, calculating the optimal solution uminThe epipolar lines corresponding to the lower left X-ray image and the right X-ray image respectively, and then the intersection point X' with the origin perpendicular to the epipolar lines respectively is calculated "1、x”2And to the intersection point x "1、x”2Performing rotation transformation to obtain the optimal projection point x'1And optimal projected point x'2Position coordinates of (2):
x’1＝T-1 1RT 1x”1(15)
x’2＝T-1 2RT 2x”2(16)。
as a specific embodiment, the process of estimating the position coordinates of the feature points in the three-dimensional space in step three (S300) includes the steps of:
step three, setting optimal projection point x'1＝(a1，b1，1)TOptimum projected point x'2＝(a2，b2，1)TOptimum projected point x'1Optimal projected point x'2With the feature point X in three-dimensional space, the camera matrix P1A camera matrix P2There is the following relationship between:
Figure BDA0002465188280000081
unfolding the above formula to obtain
Figure BDA0002465188280000082
Wherein, Pn iT(i-1, 2, 3; n-1, 2) is PnThe row(s).
The system of linear equations expressed by the equation (18) may be expressed as a homogeneous linear equation AX of 0, wherein,
Figure BDA0002465188280000083
step three, calculating a least square solution of a homogeneous equation set AX (0) by using a singular value decomposition method, wherein the position coordinate of the characteristic point X is a column of which V corresponds to the minimum characteristic value, wherein A is UDVTIs the singular value decomposition of a.
As a specific embodiment, in step four (S400), the distal femur of the three-dimensional model is spatially transformed so that the first plane is parallel to the second plane, and a quaternion q for transforming the distal femur into the projection posture is calculated1Comprises the following steps:
step four, setting the coordinates of three characteristic points on the distal femur in the three-dimensional model under a world coordinate system as X respectively1、X2、X3And the position coordinates of the three estimated characteristic points on the distal femur under the world coordinate system are X'1、X’2、X’3Vector x12、x13Is defined by a feature point X on the first plane1、X2、X3Vector of construction, vector x'12、x’13Is a feature point X 'on the second plane'1、X’2、X’3A constructed vector;
step four, calculating a normal vector n of the first plane1Normal vector n to the second plane2Respectively as follows:
Figure BDA0002465188280000091
step four and step three, calculating a vector n vertical to the normal1And the normal vector n2Unit vector n of12Sum normal vector n1And the normal vector n2The included angle of (A):
Figure BDA0002465188280000092
step four, calculating the quaternion q for transforming the distal end of the femur into the projection attitude1Comprises the following steps:
Figure BDA0002465188280000093
calculating quaternion q from proximal tibial transition to projection pose2And calculating the quaternion q of the distal femur converted to the projection attitude1The calculation method is the same, and is not described herein again.
As a specific embodiment, step five (S500) is based on the quaternion q1And quaternion q2The process of calculating the knee joint posture angle comprises the following steps:
step five, utilizing quaternion q1Representing rotation to femur coordinate system, determining final posture of femur, and using quaternion q2The rotation is represented to the tibial coordinate system and the final pose of the tibia is determined.
Let the coordinate of a point in the original femur coordinate system be xStrand of paperThe coordinates in the rotated femoral coordinate system are x'Strand of paperThe coordinate in the primary tibial coordinate system is xShinThe coordinates in the tibial coordinate system after rotation are x'ShinThen, then
Figure BDA0002465188280000094
Wherein x isStrand of paper、x’Strand of paper、xShin、x’ShinIn the above formula, the form is an imaginary four-element number;
step two, measuring the kinematics of the tibiofemoral joint, and fixing the tibia to observe the movement of the femur by taking the tibia as a reference datum;
calculating and solving a space transformation matrix M between a femur coordinate system and a tibia coordinate system after rotation by adopting a singular value decomposition method according to the position relation among a plurality of points4×4：
x’Shin＝M4×4x’Strand of paper(24)
Wherein,
Figure BDA0002465188280000101
Figure BDA0002465188280000102
Rthigh shinIs a rotation matrix between the femoral coordinate system and the tibial coordinate system, tThigh shinα, β and gamma are rotation angles around the axis of the femur coordinate system Z, Y, X respectively as translation vectors between the femur coordinate system and the tibia coordinate system;
step five, solving alpha, beta and gamma according to the formula (27), wherein the alpha, beta and gamma are the posture angles of the knee joint which are obtained:
Figure BDA0002465188280000103
the invention has the beneficial effects that: the knee joint posture measuring method based on the three-dimensional reconstruction of the characteristic points carries out the three-dimensional reconstruction of the six characteristic points in the femur and the tibia of the knee joint, avoids the influence of an iterative registration process and soft tissue artifacts, has simple and convenient process, reduces the calculated amount, improves the calculation efficiency and has higher accuracy.
The technical features of the embodiments described above may be arbitrarily combined, and for the sake of brevity, all possible combinations of the technical features in the embodiments described above are not described, but should be considered as being within the scope of the present specification as long as there is no contradiction between the combinations of the technical features.
The above-mentioned embodiments only express several embodiments of the present invention, and the description thereof is more specific and detailed, but not construed as limiting the scope of the invention. It should be noted that, for a person skilled in the art, several variations and modifications can be made without departing from the inventive concept, which falls within the scope of the present invention. Therefore, the protection scope of the present patent shall be subject to the appended claims.

Claims (7)

1. A knee joint posture measuring method based on feature point three-dimensional reconstruction is characterized by comprising the following steps:
step one, carrying out CT scanning on the knee joint, and establishing a three-dimensional model of the knee joint; calibrating an X-ray imaging system with a left X-ray source, a right X-ray source and a corresponding imaging plane to obtain a spatial position parameter of the X-ray imaging system, and establishing a virtual imaging space in three-dimensional modeling software by using the obtained spatial position parameter;
selecting six characteristic points in the three-dimensional model, and extracting position information of projection points corresponding to the six characteristic points from a left X-ray image corresponding to the left X-ray source and a right X-ray image corresponding to the right X-ray source respectively;
thirdly, forming projection rays by using the projection points in the left X-ray source and the left X-ray image and the projection points in the right X-ray source and the right X-ray image, and estimating the position coordinates of the feature points in a three-dimensional space;
step four, forming three characteristic points on the distal femur in the three-dimensional model into a first plane, forming the three characteristic points into a second plane according to the estimated position coordinates of the three characteristic points on the distal femur, carrying out spatial transformation on the distal femur of the three-dimensional model to enable the first plane and the second plane to be parallel, and calculating to obtain a quaternion q for converting the distal femur into a projection posture1；
Forming three characteristic points on the proximal tibia end in the three-dimensional model into a third plane, forming the three characteristic points into a fourth plane according to estimated position coordinates of the three characteristic points on the proximal tibia end, performing spatial transformation on the proximal tibia end of the three-dimensional model to enable the third plane and the fourth plane to be parallel, and calculating to obtain a quaternion q converted from the proximal tibia end to a projection posture2；
Step five, according to the quaternion q1And said quaternion q2And calculating to obtain a knee joint posture angle, and determining the posture of the knee joint.
2. The knee joint posture measurement method based on the feature point three-dimensional reconstruction as claimed in claim 1, wherein the process of obtaining the spatial position parameter of the X-ray imaging system in the first step comprises the following steps:
placing the left X-ray source and the right X-ray source at an included angle of 30 degrees, and performing three-dimensional calibration on the X-ray imaging system by using a Zhang-friend calibration method to obtain an internal reference matrix K of the left X-ray source and a corresponding imaging plane1External reference matrix M1A camera matrix P1And an internal reference matrix K of the right X-ray source and the corresponding imaging plane2External reference matrix M2A camera matrix P2The position relationship of the left X-ray source and the right X-ray source can be expressed by a matrix E as:
Figure FDA0002465188270000021
wherein x islIs the coordinate, X, of a certain point X in space under the camera coordinate system of the left X-ray sourcerThe coordinate of a certain point X in the space under the camera coordinate system of the right X-ray source is shown, R is a rotation matrix of the camera coordinate system of the right X-ray source and the camera coordinate system of the left X-ray source, and t is a translation vector;
the basic matrix F of the X-ray imaging system is as follows:
F＝[K2t]xK2R K1 -1(2)
wherein [ K ]2t]xIs an antisymmetric matrix of three-dimensional column vectors.
3. The knee joint posture measuring method based on the feature point three-dimensional reconstruction as claimed in claim 1, wherein the second step comprises the following steps:
selecting the medial epicondyle, the lateral epicondyle and the femoral intercondylar concavity on the distal femur and the lateral tibial condyle, the medial tibial condyle and the intercondylar eminence on the proximal tibia in the three-dimensional model as the characteristic points;
manually labeling the projection point positions of the feature points in the left X-ray image and the right X-ray image respectively;
and carrying out binarization processing on the left X-ray image and the right X-ray image after the projection point position is marked, and respectively extracting the position information of the projection point.
4. The knee joint posture measurement method based on the feature point three-dimensional reconstruction as claimed in claim 3, wherein the second step further comprises a process of calibrating the position information of the extracted projection point to obtain the position coordinates of the optimal projection point, the process comprising the steps of:
step two, assuming that the projection point extracted from the left X-ray image is X1＝(m，n，1)TThe projection point extracted from the right X-ray image is X2＝(m’，n’，1)TProjection point x1＝(m，n，1)TThe corresponding optimal projection point is x'1Projection point x2＝(m’，n’，1)TThe corresponding optimal projection point is x'2Projection point x1＝(m，n，1)TProjection point x2＝(m’，n’，1)TOptimal projected point x'1And optimal projected point x'2The distance minimization formula is satisfied:
C(x1，x2)＝d(x1、x’1)2+d(x’2，x’2)2(3)
step two, defining projection point x1＝(m，n，1)TAnd x2＝(m’，n’，1)TAnd translating it to the origin of coordinates:
Figure FDA0002465188270000031
Figure FDA0002465188270000032
step two and step three, respectively calculating a pole e in an imaging plane corresponding to the left X-ray source1＝(ex，ey，ez)TPole e in the imaging plane corresponding to the right X-ray source2＝(e’x，e’y，e’z)TPole e1And pole e2Respectively satisfy Fe10 and FTe20, pole e1And pole e2Respectively normalized to e2 x+e2 y＝1，e’x 2+e’y 21, and constructing a rotation matrix R1、R2So that it satisfies R1e1＝(1，0，ez)T、R2e1＝(1，0，e’z)TConstructed rotation matrix R1、R2Respectively as follows:
Figure FDA0002465188270000033
step two and four, the transformed basic matrix is F' ═ R2T2 -TF T1R1 TLet f be ez、f’＝e’z、a＝F’22、b＝F’23、c＝F’32、d＝F’33F' may be represented in the form:
Figure FDA0002465188270000041
parameterizing epipolar line bundles in the left X-ray image by using a parameter u, wherein the epipolar line corresponding to the left X-ray image is set as l (u), and the epipolar line l (u) can be represented as a point (0, u, 1) passing through the left X-ray imageTWith the transformed pole R in the left X-ray image1e1The cross product of (a), namely:
l(u)＝(0，u，1)T×(1，0，ez)T
＝(0，u，1)T×(1，0，f)T
＝(uf，1，-u)T(8)
let the epipolar line corresponding to the right X-ray image be l ' (u), and the epipolar line l ' (u) can be calculated by using the transformed fundamental matrix F ':
l’(u)＝F’(0，u，1)T＝(-f’(cu+d)，au+b，cu+d)T(9)
step two, projection point x1＝(m，n，1)TAnd projection point x2＝(m’，n’，1)TThe distances to the corresponding epipolar lines in the left and right X-ray images, respectively, can be expressed as:
Figure FDA0002465188270000042
Figure FDA0002465188270000043
step two, because the optimal projection point is on the epipolar line, the formula (3) is expressed as a distance function parameterized by the parameter u as:
Figure FDA0002465188270000044
the derivation of equation (12) yields:
Figure FDA0002465188270000045
the denominators are combined and the numerator is made equal to 0 to obtain
g(u)＝u((au+b)2+f'2(cu+d)2)2-(1+f2u2)2(ad-bc)(au+b)(cu+d)＝0 (14)
Step two eight, solving equation (14) to obtain 6 roots, comparing values of s (u) when u → ∞ is obtained by substituting 6 roots into equation (12), and obtaining optimal solution u (u) with s (u) taking minimum valuemin；
Step two and nine, calculating the optimal solution uminThe epipolar lines corresponding to the lower left X-ray image and the right X-ray image respectively are calculated, and the original points are perpendicular to the epipolar lines respectivelyIntersection x of lines "1、x”2And to the intersection point x "1、x”2Performing rotation transformation to obtain the optimal projection point x'1And optimal projected point x'2Position coordinates of (2):
x’1＝T-1 1RT 1x”1(15)
x’2＝T-1 2RT 2x”2(16)。
5. the knee joint posture measurement method based on the feature point three-dimensional reconstruction is characterized in that the process of estimating the position coordinates of the feature point in the three-dimensional space in the step three comprises the following steps:
step three, setting optimal projection point x'1＝(a1，b1，1)TOptimum projected point x'2＝(a2，b2，1)TOptimum projected point x'1Optimal projected point x'2With the feature point X in three-dimensional space, the camera matrix P1A camera matrix P2There is the following relationship between:
Figure FDA0002465188270000051
unfolding the above formula to obtain
Figure FDA0002465188270000052
Wherein, Pn iT(i-1, 2, 3; n-1, 2) is PnThe row(s).
The system of linear equations expressed by the equation (18) may be expressed as a homogeneous linear equation AX of 0, wherein,
Figure FDA0002465188270000061
step three, calculating the homogeneity by using a singular value decomposition methodThe least square solution of the equation set AX is 0, and the position coordinate of the characteristic point X is V corresponding to the column of the minimum characteristic value, wherein A is UDVTIs the singular value decomposition of a.
6. The knee joint posture measurement method based on feature point three-dimensional reconstruction as claimed in claim 5, wherein in step four, the distal femur of the three-dimensional model is spatially transformed so that the first plane is parallel to the second plane, and a quaternion q of the distal femur transformed to the projection posture is calculated1Comprises the following steps:
step four, setting the coordinates of three characteristic points on the distal femur in the three-dimensional model under a world coordinate system as X respectively1、X2、X3And the position coordinates of the three estimated characteristic points on the distal femur under the world coordinate system are X'1、X’2、X’3Vector x12、x13Is defined by a feature point X on the first plane1、X2、X3Vector of construction, vector x'12、x’13Is a feature point X 'on the second plane'1、X’2、X’3A constructed vector;
step four, calculating a normal vector n of the first plane1Normal vector n to the second plane2Respectively as follows:
Figure FDA0002465188270000062
step four and step three, calculating a vector n vertical to the normal1And the normal vector n2Unit vector n of12Sum normal vector n1And the normal vector n2The included angle of (A):
Figure FDA0002465188270000063
step four, calculating the quaternion q for transforming the distal end of the femur into the projection attitude1Comprises the following steps:
Figure FDA0002465188270000064
。
7. the knee joint posture measuring method based on feature point three-dimensional reconstruction of claim 6, characterized in that step five is based on the quaternion q1And said quaternion q2The process of calculating the knee joint posture angle comprises the following steps:
step five, utilizing quaternion q1Representing rotation to femur coordinate system, determining final posture of femur, and using quaternion q2Representing rotation to a tibia coordinate system, and determining the final posture of the tibia;
let the coordinate of a point in the original femur coordinate system be xStrand of paperThe coordinates in the rotated femoral coordinate system are x'Strand of paperThe coordinate in the primary tibial coordinate system is xShinThe coordinates in the tibial coordinate system after rotation are x'ShinThen, then
Figure FDA0002465188270000071
Wherein x isStrand of paper、x’Strand of paper、xShin、x’ShinIn the above formula, the form is an imaginary four-element number;
step two, measuring the kinematics of the tibiofemoral joint, and fixing the tibia to observe the movement of the femur by taking the tibia as a reference datum;
calculating and solving a space transformation matrix M between a femur coordinate system and a tibia coordinate system after rotation by adopting a singular value decomposition method according to the position relation among a plurality of points4×4:
x’Shin＝M4×4x’Strand of paper(24)
Wherein,
Figure FDA0002465188270000072
Figure FDA0002465188270000081
Rthigh shinIs a rotation matrix between the femoral coordinate system and the tibial coordinate system, tThigh shinα, β and gamma are rotation angles around the axis of the femur coordinate system Z, Y, X respectively as translation vectors between the femur coordinate system and the tibia coordinate system;
step five, solving alpha, beta and gamma according to the formula (27), wherein the alpha, beta and gamma are the knee joint attitude angles:
Figure FDA0002465188270000082
CN202010331744.7A 2020-04-24 2020-04-24 Knee joint posture measurement method based on three-dimensional reconstruction of feature points Active CN111568456B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010331744.7A CN111568456B (en) 2020-04-24 2020-04-24 Knee joint posture measurement method based on three-dimensional reconstruction of feature points

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010331744.7A CN111568456B (en) 2020-04-24 2020-04-24 Knee joint posture measurement method based on three-dimensional reconstruction of feature points

Publications (2)

Publication Number Publication Date
CN111568456A true CN111568456A (en) 2020-08-25
CN111568456B CN111568456B (en) 2023-07-14

Family

ID=72125098

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010331744.7A Active CN111568456B (en) 2020-04-24 2020-04-24 Knee joint posture measurement method based on three-dimensional reconstruction of feature points

Country Status (1)

Country Link
CN (1) CN111568456B (en)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113545956A (en) * 2021-06-01 2021-10-26 深圳大学 Design method of personalized and customized high-compliance knee joint exoskeleton
CN113689406A (en) * 2021-08-24 2021-11-23 北京长木谷医疗科技有限公司 Knee joint femoral posterior condylar point identification method and system based on motion simulation algorithm
CN117671221A (en) * 2024-02-01 2024-03-08 江苏一影医疗设备有限公司 Data correction method, device and storage medium based on knee joint limited angle image

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2008078643A1 (en) * 2006-12-25 2008-07-03 Softcube Co., Ltd. Artificial joint replacement assisting device, artificial joint replacement assisting method using same, and assisting system
CN102607526A (en) * 2012-01-03 2012-07-25 西安电子科技大学 Target posture measuring method based on binocular vision under double mediums
CN102940530A (en) * 2012-11-16 2013-02-27 昆明医科大学第一附属医院 Method for virtually building anterior cruciate ligament on femur and tibia tunnels
CN103512558A (en) * 2013-10-08 2014-01-15 北京理工大学 Binocular video pose measuring method for conical target and target patterns
CN105208247A (en) * 2015-09-02 2015-12-30 长春理工大学 Quaternion-based panoramic image stabilizing method
CN107680159A (en) * 2017-10-16 2018-02-09 西北工业大学 A kind of space non-cooperative target three-dimensional rebuilding method based on projection matrix
CN110030979A (en) * 2019-03-11 2019-07-19 北京空间飞行器总体设计部 A kind of space non-cooperative target relative pose measurement method based on sequence image
WO2019180745A1 (en) * 2018-03-21 2019-09-26 Karade Vikas Systems and methods for obtaining 3-d images from x-ray information for deformed elongate bones

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2008078643A1 (en) * 2006-12-25 2008-07-03 Softcube Co., Ltd. Artificial joint replacement assisting device, artificial joint replacement assisting method using same, and assisting system
CN102607526A (en) * 2012-01-03 2012-07-25 西安电子科技大学 Target posture measuring method based on binocular vision under double mediums
CN102940530A (en) * 2012-11-16 2013-02-27 昆明医科大学第一附属医院 Method for virtually building anterior cruciate ligament on femur and tibia tunnels
CN103512558A (en) * 2013-10-08 2014-01-15 北京理工大学 Binocular video pose measuring method for conical target and target patterns
CN105208247A (en) * 2015-09-02 2015-12-30 长春理工大学 Quaternion-based panoramic image stabilizing method
CN107680159A (en) * 2017-10-16 2018-02-09 西北工业大学 A kind of space non-cooperative target three-dimensional rebuilding method based on projection matrix
WO2019180745A1 (en) * 2018-03-21 2019-09-26 Karade Vikas Systems and methods for obtaining 3-d images from x-ray information for deformed elongate bones
CN110030979A (en) * 2019-03-11 2019-07-19 北京空间飞行器总体设计部 A kind of space non-cooperative target relative pose measurement method based on sequence image

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
AMIN AHMADI ET AL.: "3D Human Gait Reconstruction and Monitoring Using Body-Worn Inertial Sensors and Kinematic Modeling", 《IEEE SENSORS JOURNAL》 *
何巍等: "基于点云的膝关节胫骨三维CT与MRI图像配准", 《长春理工大学学报》 *

Cited By (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113545956A (en) * 2021-06-01 2021-10-26 深圳大学 Design method of personalized and customized high-compliance knee joint exoskeleton
CN113545956B (en) * 2021-06-01 2023-10-20 深圳大学 Personalized and customized high-compliance knee exoskeleton design method
CN113689406A (en) * 2021-08-24 2021-11-23 北京长木谷医疗科技有限公司 Knee joint femoral posterior condylar point identification method and system based on motion simulation algorithm
CN113689406B (en) * 2021-08-24 2022-04-08 北京长木谷医疗科技有限公司 Knee joint femoral posterior condylar point identification method and system based on motion simulation algorithm
CN117671221A (en) * 2024-02-01 2024-03-08 江苏一影医疗设备有限公司 Data correction method, device and storage medium based on knee joint limited angle image
CN117671221B (en) * 2024-02-01 2024-05-03 江苏一影医疗设备有限公司 Data correction method, device and storage medium based on knee joint limited angle image

Also Published As

Publication number Publication date
CN111568456B (en) 2023-07-14

Similar Documents

Publication Publication Date Title
US20240130792A1 (en) Orthopedic Fixation With Imagery Analysis
CN111568456B (en) Knee joint posture measurement method based on three-dimensional reconstruction of feature points
Cappello et al. Soft tissue artifact compensation in knee kinematics by double anatomical landmark calibration: performance of a novel method during selected motor tasks
US7804998B2 (en) Markerless motion capture system
Mündermann et al. The evolution of methods for the capture of human movement leading to markerless motion capture for biomechanical applications
US7929745B2 (en) Method and system for characterization of knee joint morphology
US20210012492A1 (en) Systems and methods for obtaining 3-d images from x-ray information for deformed elongate bones
CN113116523B (en) Orthopedic surgery registration device, terminal equipment and storage medium
CN117745978B (en) Simulation quality control method, equipment and medium based on human body three-dimensional reconstruction algorithm
Debarba et al. Augmented reality visualization of joint movements for physical examination and rehabilitation
Fohanno et al. Improvement of upper extremity kinematics estimation using a subject-specific forearm model implemented in a kinematic chain
Cai et al. A binocular-camera-assisted sensor-to-segment alignment method for inertial sensor-based human gait analysis
Fortin et al. A 3D visualization tool for the design and customization of spinal braces
CN111658142A (en) MR-based focus holographic navigation method and system
KR20160057024A (en) Markerless 3D Object Tracking Apparatus and Method therefor
Strydom et al. Anatomical joint measurement with application to medical robotics
CN116485850A (en) Real-time non-rigid registration method and system for surgical navigation image based on deep learning
CN114469079B (en) Body joint measurement method using LightHouse
CN116172590A (en) Knee joint movement analysis method and system
TWI441106B (en) An image based surgery training system
TWI786667B (en) Method and device for generating three-dimensional image data of human body skeletal joints
Mündermann et al. Markerless motion capture for biomechanical applications
Yin et al. AR-Based Surgical Navigation System Based on Dual Plane Calibration and 2D/3D Registration
Khameneh et al. Highly accurate automated patient-specific 3D bone pose and scale estimation using bi-planar pose-invariant patches in a CNN-Based 3D/2D registration framework
Hu et al. Robust Oral Localization Based on Multicamera Tracking of Self-Identifying Markers

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
GR01 Patent grant
GR01 Patent grant