1.一种基于特征点三维重建的膝关节姿态测量方法,其特征在于,包括以下步骤:步骤一、对膝关节进行CT扫描,建立膝关节的三维模型;对具有左X射线源、右X射线源与对应成像平面的X线成像系统进行标定,获得X线成像系统的空间位置参数,利用获得的空间位置参数在三维建模软件中建立虚拟成像空间;
步骤二、在所述三维模型中选取六个特征点,并分别在所述左X射线源对应的左X线图像和所述右X射线源对应的右X线图像中提取六个所述特征点对应的投影点的位置信息;
步骤三、利用所述左X射线源与所述左X线图像中的投影点以及所述右X射线源与所述右X线图像中的投影点组成投影射线,估计所述特征点在三维空间中的位置坐标;
步骤四、将所述三维模型中股骨远端上的三个特征点组成第一平面,根据估计得到的股骨远端上的三个特征点的位置坐标将该三个特征点组成第二平面,对所述三维模型的股骨远端进行空间变换,使所述第一平面与所述第二平面平行,计算得到股骨远端转换到投影姿态的四元数q1;
将所述三维模型中胫骨近端上的三个特征点组成第三平面,根据估计得到的胫骨近端上的三个特征点的位置坐标将该三个特征点组成第四平面,对所述三维模型的胫骨近端进行空间变换,使所述第三平面与所述第四平面平行,计算得到胫骨近端转换到投影姿态的四元数q2;
步骤五、根据所述四元数q1和所述四元数q2计算得到膝关节姿态角,确定膝关节的姿态。
2.根据权利要求1所述的基于特征点三维重建的膝关节姿态测量方法,其特征在于,步骤一中获得X线成像系统的空间位置参数过程包括以下步骤:将所述左X射线源和所述右X射线源呈30°夹角放置,利用张正友标定法对所述X线成像系统进行立体标定,获得所述左X射线源与对应的成像平面的内参矩阵K1、外参矩阵M1、相机矩阵P1和所述右X射线源与对应的成像平面的内参矩阵K2、外参矩阵M2、相机矩阵P2,所述左X射线源和所述右X射线源的位置关系可用矩阵E表示为:其中,xl为空间中某点X在所述左X射线源的相机坐标系下的坐标,xr为空间中某点X在所述右X射线源的相机坐标系下的坐标,R为所述右X射线源的相机坐标系与所述左X射线源的相机坐标系的旋转矩阵,t为平移向量;
所述X线成像系统的基本矩阵F为:
F=[K2 t]xK2R K1-1 (2)
其中,[K2 t]x为三维列向量的反对称矩阵。
3.根据权利要求1所述的基于特征点三维重建的膝关节姿态测量方法,其特征在于,步骤二包括以下步骤:选取所述三维模型中股骨远端上的股骨内上髁、股骨外上髁、股骨髁间凹以及胫骨近端上的胫骨外侧髁、胫骨内侧髁和髁间隆起作为所述特征点;
分别手动标注所述左X线图像和所述右X线图像中各个所述特征点的投影点位置;
对标注投影点位置后的所述左X线图像和所述右X线图像进行二值化处理,分别提取投影点的位置信息。
4.根据权利要求3所述的基于特征点三维重建的膝关节姿态测量方法,其特征在于,步骤二还包括对提取的投影点的位置信息进行校准,获得最优投影点的位置坐标的过程,该过程包括以下步骤:步骤二一、假设在所述左X线图像中提取的投影点为x1=(m,n,1)T,在所述右X线图像中提取的投影点为x2=(m’,n’,1)T,投影点x1=(m,n,1)T对应的最优投影点为x’1,投影点x2=(m’,n’,1)T对应的最优投影点为x’2,投影点x1=(m,n,1)T、投影点x2=(m’,n’,1)T、最优投影点x’1以及最优投影点x’2满足距离最小化公式:C(x1,x2)=d(x1、x’1)2+d(x’2,x’2)2 (3)步骤二二、定义投影点x1=(m,n,1)T和x2=(m’,n’,1)T的变换矩阵,并将其平移到坐标原点:步骤二三、分别计算所述左X射线源对应的成像平面中的极点e1=(ex,ey,ez)T和所述右X射线源对应的成像平面中的极点e2=(e’x,e’y,e’z)T,极点e1和极点e2分别满足Fe1=0和FTe2=0,对极点e1和极点e2分别进行归一化,使e2x+e2y=1,e’x2+e’y2=1,并构造旋转矩阵R1、R2,使其满足R1 e1=(1,0,ez)T、R2 e1=(1,0,e’z)T,构造的旋转矩阵R1、R2分别为:步骤二四、变换后的基本矩阵为F’=R2T2-TF T1R1T,令f=ez、f’=e’z、a=F’22、b=F’23、c=F’32、d=F’33,F’可表示为以下形式:步骤二五、用参数u参数化所述左X线图像中的对极线束,设所述左X线图像对应的对极线为l(u),对极线l(u)可以表示为其通过所述左X线图像中的某点(0,u,1)T与所述左X线图像中变换后的极点R1 e1的叉积,即:l(u)=(0,u,1)T×(1,0,ez)T
=(0,u,1)T×(1,0,f)T
=(uf,1,-u)T (8)
设所述右X线图像对应的对极线为l’(u),对极线l’(u)可以利用变换后的基本矩阵F’计算得到:l’(u)=F’(0,u,1)T=(-f’(cu+d),au+b,cu+d)T (9)步骤二六、投影点x1=(m,n,1)T和投影点x2=(m’,n’,1)T分别在所述左X线图像和所述右X线图像中与对应的对极线的距离可以表示为:步骤二七、由于最优投影点在对极线上,因此将式(3)表示为利用参数u参数化的距离函数为:对式(12)进行求导,得:
合并分母并使分子等于0,得
g(u)=u((au+b)2+f'2(cu+d)2)2-(1+f2u2)2(ad-bc)(au+b)(cu+d)=0 (14)步骤二八、解式(14)得到6个根,比较将6个根分别代入式(12)并且u→∞时s(u)的值,得到使s(u)取最小值的最优解umin;
步骤二九、计算最优解umin下左X线图像和右X线图像分别对应的对极线,然后计算原点分别垂直于对极线的交点x”1、x”2,并对交点x”1、x”2做旋转变换,得到最优投影点x’1和最优投影点x’2的位置坐标:x’1=T-11 RT1 x”1 (15)
x’2=T-12 RT2 x”2 (16)。
5.根据权利要求4所述的基于特征点三维重建的膝关节姿态测量方法,其特征在于,步骤三中估计所述特征点在三维空间中的位置坐标的过程包括以下步骤:步骤三一、设最优投影点x’1=(a1,b1,1)T,最优投影点x’2=(a2,b2,1)T,最优投影点x’1、最优投影点x’2与三维空间中的特征点X、相机矩阵P1、相机矩阵P2之间存在如下关系:展开上式,得
其中,PniT(i=1,2,3;n=1,2)为Pn的行。
式(18)表示的线性方程组可以表示为齐次线性方程AX=0,其中,步骤三二、利用奇异值分解法计算齐次方程组AX=0的最小二乘解,特征点X的位置坐标为V对应于最小特征值的列,其中A=UDVT是A的奇异值分解。
6.根据权利要求5所述的基于特征点三维重建的膝关节姿态测量方法,其特征在于,步骤四中,对所述三维模型的股骨远端进行空间变换,使所述第一平面与所述第二平面平行,计算得到股骨远端转换到投影姿态的四元数q1的过程包括以下步骤:步骤四一、设所述三维模型中股骨远端上的三个特征点在世界坐标系下的坐标分别为X1、X2、X3,估计得到的股骨远端上的三个特征点在世界坐标系下的位置坐标为X’1、X’2、X’3,向量x12、x13为第一平面上由特征点X1、X2、X3构成的向量,向量x’12、x’13为第二平面上由特征点X’1、X’2、X’3构成的向量;
步骤四二、计算第一平面的法向量n1与第二平面的法向量n2分别为:步骤四三、计算垂直于法向量n1与法向量n2的单位向量n12和法向量n1与法向量n2的夹角:步骤四四、计算股骨远端转换到投影姿态的四元数q1为:
。
7.根据权利要求6所述的基于特征点三维重建的膝关节姿态测量方法,其特征在于,步骤五根据所述四元数q1和所述四元数q2计算得到膝关节姿态角的过程包括以下步骤:步骤五一、利用四元数q1对股骨坐标系表示旋转,确定股骨的最终姿态,利用四元数q2对胫骨坐标系表示旋转,确定胫骨的最终姿态;
设一点在原股骨坐标系中的坐标为x股,在旋转后的股骨坐标系中的坐标为x’股,在原胫骨坐标系中的坐标为x胫,在旋转后的胫骨坐标系中的坐标为x’胫,则其中,x股、x’股、x胫、x’胫在上式中为虚四元数形式;
步骤五二、胫股关节运动学测量,以胫骨为参考基准,固定胫骨观测股骨的运动;
根据多个点之间的位置关系,采用奇异值分解法计算求取旋转后股骨坐标系与胫骨坐标系之间的空间变换矩阵M4×4:x’胫=M4×4x’股 (24)
其中,
R股胫为股骨坐标系与胫骨坐标系之间的旋转矩阵,t股胫为股骨坐标系与胫骨坐标系之间的平移向量,α、β、γ分别为绕股骨坐标系Z、Y、X轴的旋转角度;
步骤五三、根据式(27)求解α、β、γ,α、β、γ即为膝关节姿态角: