CT影像中骨的运动和接近参数的标定和量化方法与流程

本发明属于临床医学、运动医学和计算机图形学领域，涉及一种ct影像中骨的运动和接近参数的标定和量化方法。

技术背景

体育运动和日常劳作经常会导致关节中韧带的损伤，使韧带不能有效地限制骨的活动范围，进而造成关节运动不稳；若不能及时诊断和治疗会导致骨的损伤，并诱发关节炎。临床上可以通过抽屉试验、沃森试验、关节镜等手段来诊断关节不稳，也可以通过x光、ct和mr等影像学方法诊断。影像学方法具有可量化测量的优点，因此在关节不稳的诊断中具有重要作用。现有的影像学方法大多针对关节在特定姿态下的静态成像，通过测量特定骨之间距离和夹角的异常来诊断关节不稳。这种基于静态成像的方法对于关节不稳的诊断是有局限性的。动态型关节不稳(dynamicinstability)的影像学特征发生在关节的运动过程中，在静态成像方法的几个预设关节姿态上并无影像学异常，因此这类关节不稳常常被漏诊。

四维ct的出现为解决这个问题带来了新的机遇。与传统ct相比，四维ct成像速度更快、单次成像辐射剂量更低。近年来四维ct技术趋于成熟，以320排动态容积成像技术为代表的新一代ct可以实现无相差、大扫描范围的四维ct成像，生成扫描直径16公分，每秒3-6帧的运动影像视频。与其他可以动态成像的医学影像设备相比，四维ct的突出优点是能够对骨性成分进行快速清晰的体成像，进而通过计算机图形学方法可从中准确地获取骨的三维形状和运动信息，比如精确测量骨的位移和骨间距等参数，这是诊断关节动态不稳的关键依据。

利用四维ct形成的“视频”影像来量化骨的运动信息和接触信息，并用于诊断关节动态不稳将会是一种重要辅助诊断手段。目前基于四维ct的关节不稳研究，主要从影像中采集运动和接近两方面的量化参数。其中运动参数包括：骨的三维平移和旋转信息，骨的接近面中心点在骨表面的轨迹信息；接近参数包括：骨接近区域、及其与相邻骨的最小间距和平均间距。然而在这些量化参数中，只有接近区域与相邻骨的最小间距和平均间距分别被用于诊断关节不稳，其他参数只被用来研究关节运动的一般特性，并没有与关节不稳的诊断建立直接联系。这些未被使用的参数与关节不稳之间是否存在关联是需要进一步研究的；但是目前这些参数用于关节不稳的诊断是有困难的，主要原因是：这些参数是没有经过标定的，难以进行不同患者之间的量化对比。

为了解决关节骨运动参数标定问题，需要一种可以将运动和接近参数进行归一化和精确对比测算的方法，使得骨的三维运动信息、接近面和接近面中心点等参数可以在多患者的影像分析结果之间进行对比分析。

技术实现要素：

本发明的目的是，为解决骨的位置和姿态信息参数标定问题，提供一种ct影像中骨的运动和接近参数的标定和量化方法，该方法通过空间坐标变换方法使得多样本间骨具有一致性的局部坐标系，在此基础上精确计算骨的运动参数和接近参数，并进行对比，可辅助医生对关节动态不稳的诊断，减少误诊漏诊的发生。

本发明的技术解决方案如下：

一种ct影像中骨的运动参数标定和量化方法，包括以下步骤：

首先，采集样本ct影像，对该影像样本进行骨的分割和跟踪，从不同样本ct影像中提取出骨的三维模型和运动参数作为数据源；然后，选取其中一个样本ct影像中的骨作为基准骨，在该基准骨上创建局部坐标系；再将其他样本的骨与基准骨对齐，使得所有样本的骨都具有与基准骨一致性的局部坐标系；最后，基于骨的局部坐标系，对单个骨进行定位定姿。

一种基于ct影像骨的运动参数标定方法，具体包括以下步骤：

步骤1、定义一个起基准作用的骨，即基准骨b^ref，基于骨的三维模型建立局部坐标系，记为其中和分别代表基准坐标系c^ref的坐标原点及三个轴向，其所属的图像为基准图像；

步骤2、对于非基准图像，设为样本k的第t帧图像，通过骨的分割和跟踪获得其相应的骨的分割图像b^k|t，通过求解以下公式求解将b^k|t与基准骨b^ref进行骨的对齐的问题：

其中，argminθ′表示求解使得评估函数值最小的θ′值；函数是图像匹配程度的评估函数，常用的包括互信息量和图像残差均方值等；函数ψ′表示将骨的分割图像b^k|t分别进行变换矩阵为mr的三维自旋变换和变换矩阵为mt的三维平移变换，然后进行s倍的缩放；θ′表示mt和mr这两个变换矩阵中的未知数和缩放系数s；

步骤3、通过求解出的参数θ′和以下公式估算样本k的第t帧图像中骨的局部坐标系其中和分别代表c^k|t的坐标原点及三个轴向：

其中，和分别代表步骤2中三维自旋变换和三维平移变换的逆变换；

步骤4、通基于骨的局部坐标系，对单个骨进行定位定姿。

进一步地，所述骨的分割和跟踪采用手动、半自动或者全自动的分割跟踪方法

所述公式一的求解方法为：首先通过通用的医学图像处理软件来实现基准骨和非基准骨的配准；然后由非基准骨的初始平移、旋转和缩放信息以及配准后的新的平移、旋转和缩放信息推算出三维自旋变换mr和三维平移变换mt以及缩放系数s。

由于两个骨初始的位置大小以及旋转角度不同，因此配准过程包含了三维自旋变换mr和三维平移变换mt以及缩放s。由于设定了基准骨的局部坐标系其中和分别代表基准坐标系c^ref的坐标原点及三个轴向，因此通过mr和mt的逆变换可推算出与其配准的非基准骨的局部坐标系具体的实施步骤是如下：

选用一款通用医学图像处理软件对基准骨图像和非基准骨图像进行模板配准，通过模板图像的配准位移和旋转角度，可以推算出基准骨和非基准骨在配准过程中的三维自旋变换mr和三维平移变换mt以及缩放s。例如采用amira软件实现基准骨和非基准骨的配准，跟踪计算的步骤如下：

1.用amira打开基准骨图像；

2.用amira打开非基准骨图像；

3.添加仿射配准(affineregistration)节点，将其model输入端与基准骨图像节点相连，将其reference输入端与非基准骨图像节点相连；

4.在affineregistration节点的属性面板上设置变换类型为rigid(刚性配准)和iso-(伸缩配准)；

5.在非基准骨图像节点的属性面板中点击transformeditor按钮，以激活模板图像的交互式空间变换功能，将非基准骨图像进行适当移动和旋转，使其与基准骨图像大致对齐；

6.在affineregistration节点的属性面板上点击register按钮，进行模板配准；

7.在三维视图中观察配准后的非基准骨图像是否与基准骨图像精确对齐；如果是，读取并记录非基准骨图像经过配准后新的平移、旋转信息以及缩放信息。

由于非基准骨的初始平移、旋转和缩放信息可以在未变换前读出，并根据配准后的新的平移、旋转和缩放信息推算出三维自旋变换mr和三维平移变换mt以及缩放系数s。又由于基准骨的局部坐标系是已知的，因此可以通过三维自旋变换mr和三维平移变换mt的逆变换推出非基准骨的局部坐标系即解出了公式二。

本发明还公开了一种基于骨的运动参数标定的踝关节的运动姿态获取方法，根据上述方法得到胫骨下端的局部坐标系和距骨的局部坐标系；以胫骨下端的局部坐标系为基准，将其视作踝关节的基准坐标系并计算距骨的局部坐标系与c^std的相对位移和偏转角度，计算公式如下：

其中，djoint是距骨与胫骨下端局部坐标系中心的归一化距离，该值反映了踝关节运动过程中距骨的相对平移量；函数dist()用于计算两点之间的欧式距离；为踝关节的基准坐标系的原点位置；s为公式一中的缩放倍数；θjoint，ωjoint分别为距骨局部坐标系x，y，z三个轴与踝关节的基准坐标系对应轴的夹角，如图4所示，其值是通过向量点积的反余弦计算的；指距骨x轴与踝关节的基准坐标系x轴的向量点积，|·|是指相应向量的模；

定义三元组为踝关节的运动姿态，它反映了以胫骨下端为基准时距骨的偏转量，体现了两骨的相对姿态。

本发明还公开了一种ct影像中骨的接近参数标定和量化方法，根据上述方法得到骨的局部坐标系；设vi为骨的三维模型上的第i个顶点，其在骨的局部坐标系中的坐标为采用以下步骤将骨的局部坐标系上的点映射到二维骨形图上：

首先，通过c^k|t和确定vi在二维骨形图中的经纬度坐标(lngi,lati)：

其中，函数m用于将局部坐标系c^k|t下的点映射二维骨图平面上，λ代表映射的方法，可采用任何一种现有的地图投影方法。本发明选择墨卡托投影法(mercatorcylindricalprojection)作为映射方法λ，能够保持方向和角度的正确性。

然后，将样本k的第1帧图像中骨的局部坐标系c^k|1的坐标原点作为计算高程信息的基点，计算vi在二维骨形图中的高程信息，计算公式为：

其中，hlng,lat是二维骨形图上坐标为(lng,lat)的点对应高程信息；s^k|1是样本k的第1帧图像对应的骨的三维模型，函数dist()用于计算两点之间的欧式距离，函数ray(v1,v2)用于计算以v1为起点指向v2的三维射线；函数intersect()用于计算三维射线与三维模型之间的交点；对于不封闭的骨，设置高程阈值，当intersect函数无法测得交点时，将hlng,lat设置为该高程阈值；在后续的可视化或对比分析中忽略该高程阈值所在的骨图位置；s为公式一中的缩放倍数。经度lng的取值范围是[-180,180]，纬度lat的取值范围是[-90,90]；对于不同的高程信息hlng,lat进行着色处理，获得二维骨形图；

最后，以接近面s^talus顶点集合所组成的三角形网格区域的面积来估算接近面面积；

将接近面中心点和最近接触点映射到骨形图中分析其轨迹特征；最近接触点的映射方法如公式四所示；接近面中心点的计算则是通过将接近面顶点集合中的轮廓顶点映射到骨形图上，以其映射后的二维轮廓的形状重心作为接近面中心点。

本发明的原理为：

1.基于局部坐标系的运动参数标定方法。骨是一类具有相似形状的三维物体，它们外形呈现不规则的起伏，而且起伏比较平缓，很少有可以精确定位的尖角或特征位置。因此，这类物体不能像长方体那样可以简单地根据尖角(特征位置)来定位和定姿。运动参数标定的关键就是找到能对其定位定姿的方法。

假设骨是具有相似三维外形的实心物体，它的空间属性包含四个方面：位置、姿态、尺度和外形。其中外形是指与位置、姿态、尺度无关的分布属性，在这四个属性中位置和姿态可以用骨的局部坐标系来表示，c＝{po，ax，ay，az}；其中po指骨的局坐标原点在骨的三维模型中的位置；ax，ay，az指骨的局部坐标系的三个正交轴在骨的三维模型中的朝向。

由于骨表面一般没有可以精确定位的尖角，因此无法根据其外形定义一种普适于每个样本的局部坐标系原点及轴向，即难以找到一种绝对的独立于骨样本差异性之外的局部坐标系定义方法。但可以采用相对的策略来定义，即对某个样本的骨指定其局部坐标系，然后利用形状相似性将其它样本的骨与该骨对齐，以此来定义每个样本的骨的局部坐标系。定义一个起基准作用的骨，即基准骨b^ref，其局部坐标系为c^ref，其所属的图像为基准图像i^ref；

对于非基准图像而言，虽然其骨的局部坐标系是未知的，但知道其理想的局部坐标系应该与c^ref一致，即当把这些骨图像进行变换，使其对应的骨的局部坐标系与c^ref的原点重合、轴向重合且尺度归一化，则此时不同样本的骨与基准骨之间应该实现最大程度的重叠或对齐，即只有外形的样本差异存在。只有这样的局部坐标系才是与c^ref一致的，才可起到界定其位置和姿态的作用。该思路可以用公式(1)表示。

argminθφ[ψ(c^k|t,s^k|t,i^k|t,c^ref),i^ref](1)

其中，argminθ表示求解使得比较函数值φ最小的变量值θ；φ是基于图像的骨的差异性的比较函数，ψ是图像变换方法，将样本k的第t帧图像进行倍率为s^k|t的缩放，并进行适当平移和旋转变换，使其对应的骨的局部坐标系c^k|t与基准骨的局部坐标系c^ref重叠；θ为c^k|t中的未知参数；i^k|t为样本k的第帧图像对应的骨的位置参数，i^ref为基准骨的位置参数，用于平移变换；

直接求解公式(1)是困难的，因为φ一般是非线性的、非凸的，因此其优化精度往往取决于初始值的准确性；由于拍摄条件和样本差异性的存在，直接估算c^k|t的初始值是难以准确的。可以通过手工或半自动分割获得b^k|t，即骨的分割图像。由于将b^k|t与基准骨b^ref进行精确的骨的对齐是可行的，因此可以把公式(1)中的优化问题，转换为b^k|t与基准骨b^ref这两个分割图像之间的骨的对齐问题，如公式(2)所示。

其中，argminθ′表示求解使得评估函数值最小的θ′值；函数是图像匹配程度的评估函数，常用的包括互信息量和图像残差均方值等；函数ψ′是指将骨的分割图像b^k|t分别进行变换矩阵为mr的三维自旋变换和mt的三维平移变换，然后进行s倍的缩放；θ′是指mt、mr这两个变换矩阵中的未知数和缩放系数s。

通过公式(2)的求解可以估算出样本k的第t帧图像中骨的局部坐标系如公式(3)。

其中，和分别代表公式(2)中三维自旋变换和三维平移变换的逆变换。

通过上述方法可以获得任意样本图像中骨的局部坐标系c^k|t。这个局部坐标系可以对单个骨进行定位定姿，但对于关节不稳的分析而言，更重要的是距骨和胫骨下端的相对姿态，而非单个骨的姿态。如图3所示。

在踝关节的运动中，胫骨下端一般是静止的，而距骨则围绕胫骨下端旋转。因此本发明以胫骨下端的局部坐标系为基准，将其视作踝关节的基准坐标系并计算距骨的局部坐标系与c^std的相对位移和偏转角度，以此来定义踝关节的运动姿态，如公式(4)。

其中，djoint是距骨与胫骨下端局部坐标系中心的归一化距离，该值反映了踝关节运动过程中距骨的相对平移量；函数dist()用于计算两点之间的欧式距离；为基准坐标系的原点位置；s为公式(2)中的缩放倍数。引入缩放倍数s是为了消除骨大小(尺寸)不同这一样本差异；θjoint，ωjoint分别为距骨局部坐标系x，y，z三个轴与踝关节的基准坐标系对应轴的夹角，如图4所示，其值是通过向量点积的反余弦计算的；指距骨x轴与基准坐标系x轴的向量点积，|·|是指相应向量的模。

本发明以三元组定义踝关节的姿态。它反映了以胫骨下端为基准时距骨的偏转量，体现了两骨的相对姿态。由公式(2)和(3)可知，各骨的局部坐标系在不同样本图像中是具有一致性的，因此以该三元组定义的踝关节姿态可以解决不同样本骨的坐标系不同导致运动参数不能对比的问题。

2.基于二维骨形图的接近参数标定方法。本发明提出二维骨图法，将骨的三维表面展开为具有经纬度坐标的二维图像。基于二维骨图，可对骨表面上的任意点进行定位，并对从属于该点的接近参数进行对比分析和可视化。二维骨图提供了骨表面的参数化方法，使得表面上的各项参数可定位、可对比，解决了接近面参数无法在三维骨表面精确标定问题。

把骨看作是类似地球的封闭三维形体，仿照地理测绘学中的地图映射法，把骨三维表面上的每一个点映射到二维骨图上来。这样，二维骨图上的每一个点就对应了三维骨表面上的一个点。如果还在二维骨图上标记该点在三维骨表面上对应点的高程信息，就可以根据二维骨图复原骨的三维模型。因此二维骨图与三维骨表面相比，虽然降低了维度，但并没有丢失几何信息。

设vi为骨的三维网格模型上的第i个顶点，其在骨的局部坐标系中的坐标为二维骨图计算的关键是通过c^k|t和确定vi在二维展开图像中的经纬度坐标(lngi,lati)，如公式(5)。

其中函数m用于将坐标系c^k|t下的点映射二维骨图平面上(m是一个类似地图映射函数的函数，能够将三维表面点映射到二维坐标系中)，λ代表映射的方法。选择最常见的墨卡托投影法(mercatorcylindricalprojection)，该方法的特点是能够保持方向和角度的正确，常见的世界地图就是采用就是这种映射方法。

为了将骨的表面形状进行量化对比，本发明对不同样本中的距骨创建了其对应的二维骨形图，即像包含了地形信息的地图一样，地图上的每个点既有经纬度坐标还有高程信息，以此来对地理特征进行定位和描述。其计算方法如公式(6)。

其中hlng,lat是二维骨图上坐标为(lng,lat)的点对应高程信息。由于在骨的表面上没有类似海平面的结构，因此将计算海拔的基点设置在骨的局部坐标系c^k|1的坐标原点s^k|1是骨的三维模型，函数dist()用于计算两点之间的欧式距离，函数ray(v1,v2)用于计算以v1为起点指向v2的三维射线。函数intersect()用于计算三维射线与三维模型之间的交点；本发明选择样本的第1帧对应的骨的局部坐标系的坐标原点和骨的三维模型s^k|1进行计算，这是因为骨在运动中近似不发生形变，因此各帧的骨形图应该是相同的。s为公式(2)中的缩放倍数。经度lng和纬度lat的取值范围根据选用地图映射方法进行合理取值，例如经度lng的取值范围是[-180,180]，也可以是[0,360]，纬度lat的取值范围是[-90,90]。对于不同的高程信息hlng,lat进行着色处理，可以获得如图5所示的二维图形图。

本发明以接近面s^talus顶点集合所组成的三角形网格区域的面积来估算接近面面积。对于不封闭的骨，比如胫骨下端在踝关节成像中是局部成像的，因此其三维模型是开口的。可以设置高程阈值，当intersect函数无法测得交点时，将hlng,lat设置为该高程阈值。后续的可视化或对比分析中忽略该高程阈值所在的骨图位置即可。同时将接近面中心点和最近接触点映射到骨形图中分析其轨迹特征。最近接触点的映射方法如公式(5)所示。接近面中心点的计算则是通过将接近面顶点集合中的轮廓顶点映射到骨形图上，以其映射后的二维轮廓的形状重心作为中心点。

本发明提出了一种ct影像中骨的运动和接近参数的标定和量化方法，包括：首先采集患者关节运动过程的ct影像，对该影像样本进行骨的分割和跟踪，提取骨的三维形状和运动信息；然后对多样本骨的位置和姿态进行归一化处理，即建立具有一致性的局部坐标系，在此基础上精确计算出多样本骨不同帧时的运动参数信息，可用于对比；同时在已有的一致性局部坐标系，将三维骨表面展开成二维图形，仿照地球表面展开成地图的方法，形成经纬坐标，可精确标定三维骨表面的接近信息，同样可用于对比。

有益效果：

本发明提出的ct影像中骨的运动和接近参数的标定和量化方法，可有效的解决诸如运动轨迹、接近面和中心点的位移等参数虽然可以提取但无法直接进行对比分析的问题，为ct影像分析在关节不稳诊断中的深入应用，辅助医生诊断，减少误诊漏诊的发生打下了基础。本发明通过手动分割和半自动跟踪方法以尽量保证骨的分割和跟踪精度的基础上，进行骨的运动参数和接近参数标定，准确性高。

附图说明

图1是运动参数标定的必要性示意图；其中图1(a)显示的是接近区域中心点在4个关节姿态之间的运动轨迹；图1(b)显示的是两个患者的接近区域中心点沿x坐标轴的平移曲线；图中trapezium表示大多角骨，initialkeypinchposition表示初始感到疼痛的位置，adduction表示内收，palmarabduction表示掌外展；opposition表示对立位置。由于两个坐标系不一致，这两条曲线不能直接对比。

图2是接近参数标定的必要性示意图；其中图2(a)显示的是接近区域在骨表面不同区域的分布情况，图中深色区域代表接近区域，虚线构成的网格是对骨表面的分区，从图中可见接近区域在不同关节姿态时的所处位置不同，图中的数字代表分布在该分区的接近区域占总接近区域的百分比；图中trapezium表示大多角骨，initialkeypinchposition表示初始感到疼痛的位置，dorsal表示背侧的，ulnar表示尺骨的，volar表示手掌的，radial表示桡骨，adduction表示内收；图2(b)显示的是三维物体表面的不可定位性；其中立方体表面任意点是可以精确定位的，因为立方体的任意平面都有4个可精确定位的顶点，这4个点又可以确定一个坐标系来对表面任意点定位。但球体表面上的点就无法精确定位，因为球体圆滑的表面上没有任何可定位的基准点。与球体类似，骨的表面上也缺乏可以精确定位的基准点。

图3是相对姿态对关节不稳分析的重要性示意图，左图是平放的踝关节，右图是整体抬起小腿的踝关节。虽然距骨旋转了，但距骨与胫骨下端的相对姿态没有发生改变，踝关节并没有运动。

图4是踝关节姿态参数示意图。

图5是二维骨形图示例，其中图5(a)为三维距骨模型，图5(b)为对应的骨形图。骨形图以世界地图的方式呈现，坐标为经纬度。颜色越暖的地方代表高程越大，颜色越冷的地方代表高程越小。

图6是距骨沿基准坐标系x轴的平移曲线图。

图7是距骨沿基准坐标系y轴的平移曲线图。

图8是距骨沿基准坐标系z轴的平移曲线图。

图9是最近骨间距在不同姿态下的取值图。

图10是距骨接近面的平均骨间距在不同姿态下的取值图。

图11是距骨接近面面积在不同姿态下的取值图。

图12是运动过程中的关节姿态曲线图。

图13是双侧距骨的二维骨形图；图13(a)和(b)分别为左侧和右侧距骨的二维骨形图。

具体实施方式

下面结合图示附图和具体实验过程对本发明做进一步说明。需要说明的是，本实施例中所采用的骨的分割和跟踪方法、医学图像处理软件、映射的方法不构成对本发明保护范围的限制，本发明步骤2不局限于特定的骨的分割和跟踪方法，公式一的求解不局限于特定的医学图像处理软件、公式四中的地图映射方法也不局限于特定的地图映射方法，任何通过对本发明中技术特征进行等同替换得到的技术方案均属于本发明的保护范围。

本发明公开了一种ct影像中骨的运动参数标定和量化方法，具体包括以下步骤：

步骤1、定义一个起基准作用的骨，即基准骨b^ref，基于骨的三维模型建立局部坐标系，记为其中和分别代表基准坐标系c^ref的坐标原点及三个轴向，其所属的图像为基准图像；

步骤2、对于非基准图像，设为样本k的第t帧图像，通过骨的分割和跟踪获得其相应的骨的分割图像b^k|t，通过求解以下公式求解将b^k|t与基准骨b^ref进行骨的对齐的问题：

其中，argminθ′表示求解使得评估函数值最小的θ′值；函数是图像匹配程度的评估函数，常用的包括互信息量和图像残差均方值等；函数ψ′表示将骨的分割图像b^k|t分别进行变换矩阵为mr的三维自旋变换和变换矩阵为mt的三维平移变换，然后进行s倍的缩放；θ′表示mt和mr这两个变换矩阵中的未知数和缩放系数s；

步骤3、通过求解出的参数θ′和以下公式估算样本k的第t帧图像中骨的局部坐标系其中和分别代表c^k|t的坐标原点及三个轴向：

其中，和分别代表步骤2中三维自旋变换和三维平移变换的逆变换；

步骤4、通基于骨的局部坐标系，对单个骨进行定位定姿。

所述公式一的求解方法为：首先通过医学图像处理软件来实现基准骨和非基准骨的配准；然后由非基准骨的初始平移、旋转和缩放信息以及配准后的新的平移、旋转和缩放信息推算出三维自旋变换mr和三维平移变换mt以及缩放系数s。

采用一款通用医学图像处理软件对基准骨图像和非基准骨图像进行模板配准，通过模板图像的配准位移和旋转角度，可以推算出基准骨和非基准骨在配准过程中的三维自旋变换mr和三维平移变换mt以及缩放s。

本发明还公开了一种基于骨的运动参数标定的踝关节的运动姿态获取方法，根据上述方法得到胫骨下端的局部坐标系和距骨的局部坐标系；以胫骨下端的局部坐标系为基准，将其视作踝关节的基准坐标系并计算距骨的局部坐标系与c^std的相对位移和偏转角度，计算公式如下：

其中，djoint是距骨与胫骨下端局部坐标系中心的归一化距离，该值反映了踝关节运动过程中距骨的相对平移量；函数dist()用于计算两点之间的欧式距离；为踝关节的基准坐标系的原点位置；s为公式一中的缩放倍数；θjoint，ωjoint分别为距骨局部坐标系x，y，z三个轴与踝关节的基准坐标系对应轴的夹角，如图4所示，其值是通过向量点积的反余弦计算的；指距骨x轴与踝关节的基准坐标系x轴的向量点积，||是指相应向量的模；

定义三元组为踝关节的运动姿态，它反映了以胫骨下端为基准时距骨的偏转量，体现了两骨的相对姿态。

本发明还公开了一种基于ct影像骨的接近参数标定方法，根据上述方法得到骨的局部坐标系；设vi为骨的三维模型上的第i个顶点，其在骨的局部坐标系中的坐标为采用以下步骤将骨的局部坐标系上的点映射到二维骨形图上：

首先，通过c^k|t和确定vi在二维骨形图中的经纬度坐标(lngi,lati)：

其中，函数m用于将局部坐标系c^k|t下的点映射二维骨图平面上，λ代表映射的方法。本发明选择墨卡托投影法(mercatorcylindricalprojection)作为映射方法λ，能够保持方向和角度的正确性。

然后，将样本k的第1帧图像中骨的局部坐标系c^k|1的坐标原点作为计算高程信息的基点，计算vi在二维骨形图中的高程信息，计算公式为：

其中，hlng,lat是二维骨形图上坐标为(lng,lat)的点对应高程信息；s^k|1是样本k的第1帧图像对应的骨的三维模型，函数dist()用于计算两点之间的欧式距离，函数ray(v1,v2)用于计算以v1为起点指向v2的三维射线；函数intersect()用于计算三维射线与三维模型之间的交点；对于不封闭的骨，设置高程阈值，当intersect函数无法测得交点时，将hlng,lat设置为该高程阈值；在后续的可视化或对比分析中忽略该高程阈值所在的骨图位置；s为公式一中的缩放倍数。经度lng和纬度lat的取值范围根据选用地图映射方法进行合理取值，本实施例中经度lng的取值范围是[-180,180]，纬度lat的取值范围是[-90,90]；对于不同的高程信息hlng,lat进行着色处理，获得二维骨形图；

最后，以接近面s^talus顶点集合所组成的三角形网格区域的面积来估算接近面面积；

将接近面中心点和最近接触点映射到骨形图中分析其轨迹特征；最近接触点的映射方法如公式四所示；接近面中心点的计算则是通过将接近面顶点集合中的轮廓顶点映射到骨形图上，以其映射后的二维轮廓的形状重心作为接近面中心点。

以下是针对踝关节不稳的多样本对比实验。需要说明的是，本发明虽然以四维ct为实验背景，但是本发明中的运动参数和接近参数标定方法能应用于所有ct影像。

1、实验数据

本发明选择6位志愿者的四维ct影像作为实验数据。志愿者在医师的指导下，分别将左右踝关节从跖屈极限位运动到背屈极限位。四维ct对该过程连续成像获得9-14副三维图像。其中有5位志愿者拍摄了双侧脚踝的运动影像，1位志愿者拍摄了单侧脚踝的运动影像。每位志愿者都有一侧脚踝是受伤的。表4-1显示了这6位志愿者及其四维ct影像的基本情况。本发明以志愿者2的右踝距骨和胫骨下端作为基准骨来标定其他志愿者的运动和接近面参数。

表1志愿者及其四维ct影像基本情况表

2、实验方法

采用手动分割和半自动跟踪的方法，提取距骨和胫骨下端的形状和运动参数；然后采用本发明中提出的基于ct影像的骨的运动参数和接近参数标定方法对样本参数进行标定，并进行各参数之间的量化对比。并讨论各参数与踝关节不稳的关系。

3、实验结果

在关节运动过程中，胫骨的偏移不是很明显，主要考察距骨的位置和姿态信息，通过实验得到可对比的多样本运动信息和接近信息变化趋势。本发明以三元组定义踝关节的姿态，它与其他参数构成超过3维的高维数据，因此无法通过常规的二维和三维曲线来展现不同姿态下的运动和接近参数。为了解决这一问题，采用平行坐标图(parallelcoordinatesplot)。该图可以展现高维数据内部的关系。以下是样本2经过归一化处理后可用于不同样本间对比的健康一侧脚踝和受伤一侧脚踝的运动参数和接近参数的标定结果。

(1)距骨沿基准坐标系x轴的平移曲线图如图6所示。图6中实线代表左脚的参数，虚线代表右脚的参数。第1至第3列为踝关节姿态的三元组在不同帧的样本值，单位是度。第4列为骨沿踝关节基准坐标系x轴的平移值，单位是毫米。从图中可以看出扭伤的右脚在x方向的平移相比左脚而言幅度要小，只在11-12毫米之间；而健康的左脚在x方向的平移幅度超过了13毫米。

(2)距骨沿基准坐标系y轴的平移曲线图如图7所示。图7中实线代表左脚的参数，虚线代表右脚的参数。与图6相似，不过第4列为骨沿基准坐标系y轴的平移值，单位是毫米。从图中可以看出可见扭伤的右脚在y方向的平移相比左脚而言幅度要较小，只在[0.14,1.8]毫米之间；而健康的左脚在y方向的平移幅度在[-0.4,2.2]毫米之间。

(3)距骨沿基准坐标系z轴的平移曲线图如图8所示。图8中实线代表左脚的参数，虚线代表右脚的参数。与图6相似，不过第4列为骨沿基准坐标系z轴的平移值，单位是毫米。从图中可以看出可见扭伤的右脚在z方向的平移相比左脚而言幅度要较小，只在[-29,-24]毫米之间；而健康的左脚在该z向的平移幅度在[-30,-22]毫米之间。

(4)最近骨间距在不同姿态下的取值图如图9所示。图9显示了不同姿态下最近骨间距的取值。与图6相似，不过第4列为最近骨间距，单位是毫米。从图中可以看出扭伤的右脚整体上的最近骨间距要小于健康脚。从曲线的走向来看，在前三列走向相同的曲线，在第三与第四列之间呈现不同的走势。右脚的曲线呈现下降趋势，而左脚的曲线呈现上升趋势。健康脚的最近骨间距可以达到2.475毫米，而扭伤的与右脚则只有不足1.8毫米。

(5)距骨接近面的平均骨间距在不同姿态下的取值如图10所示。图10显示了不同姿态下接近面的平均骨间距的取值。与图6相似，不过第4列为接近面的平均骨间距，单位是毫米。从图中可以看出扭伤的右脚整体上的接近面的平均骨间距方面也要小于健康脚。从曲线的走向来看，在前三列走向相同的曲线，在第三与第四列之间呈现不同的走势。右脚的曲线呈现下降趋势，而左脚的曲线呈现上升趋势。健康脚接近面的平均骨间距可以达到4毫米，而扭伤的右脚则只有不足3.5毫米。这与图9表现出的规律性是一样的。

(6)距骨接近面面积在不同姿态下的取值如图11所示。图11显示了不同姿态下接近面的面积的取值。与图6相似，不过第4列为接近面的面积，单位是平方毫米。从图中可以看出扭伤的右脚整体上的接近面的平均骨间距要大于健康脚。从曲线的走向来看，在前三列走向相同的曲线，在第三与第四列之间呈现相似的走势，都呈交叉状，即一部分曲线向下延伸，另一部分曲线向上延伸。这说明两侧脚的接近面积的变换率都经历了较快幅度的变化。健康脚接近面的面积最大为1320平方毫米，而扭伤的与右脚则可以达到1455平方毫米。这与图10表现出的规律性是一样的，因为右脚的平均骨间距小，因此距骨与胫骨下端更靠近，接近面的面积也自然更大。

(7)运动过程中的关节姿态曲线图如图12所示。图12中坐标系的三个轴分别为姿态参数三元组虚线代表左脚的参数，实线代表右脚的参数。虚线在三个方向都长于实线，尤其在ωjoint方向更加明显。实线只能下降到20多度，这说明右脚无法到达背屈极限位。可见受伤的右脚的姿态变化幅度是小于健康的左脚的。

(8)双侧距骨的二维骨形图如图13所示。在图13中左右两图分别是左侧和右侧距骨的骨形图。骨形图以世界地图的方式呈现，坐标为经纬度。颜色越暖的地方代表高程越大，颜色越冷的地方代表高程越小。