基于双平面X光的半自动配准方法和关节三维运动求解算法

基于双平面x光的半自动配准方法和关节三维运动求解算法
技术领域
1.本发明涉及三维运动捕捉分析领域，具体涉及一种基于双平面x光的半自动配准方法和基于双平面x光半自动配准数据的关节三维运动求解算法。


背景技术：

2.关节运动分析是基于运动学原理并结合解剖学和生理学及运动分析技术对人体关节运动进行定量的描述；分析关节运动规律有着重要的意义，在工程方面，对于关节植入物的设计和康复辅具的开发可提供重要的生物力学基础和依据；在临床方面，对于关节疾病的辅助诊断和术后康复的评估也具有重要意义。
3.人体各部分骨骼大小不一，形态各异，如足部骨骼、手部骨骼等，而且这些骨骼构成的关节的运动表现为空间六自由度的变化，这就对研究人体关节运动提出了严峻的挑战。
4.早期关节运动研究多采用离体试验，手动驱动离体关节样本或采用机器人驱动器驱动样本模拟特定动作，完成对相关骨骼运动的捕捉。显然，离体试验具有很大的局限性，基于运动模拟装置的试验，无法再现人体关节的真实运动模式，无法真实还原在体的关节运动状况；现代运动分析多采用远红外光运动捕捉系统(optical motion capture，omc)，omc的原理是通过在皮肤上贴标记点，使用远红外光捕捉标记点的运动，通过标记点的运动数据和人体多节段模型来解析骨骼的空间运动特征；但是，此方法也有明显的弊端：(1)这类研究存在很多前提假设和对研究问题的简化，例如将人体关节划分为几个刚性节段。(2)由于皮肤滑移的影响，皮肤表面的标记点和皮下骨骼有相对滑移运动；因此，对于骨骼较小、关节较多的人体足部或手部关节运动研究其误差更为明显；此外，有研究采用局部麻醉辅助办法在目标骨骼嵌入骨钉，且在体表设置与骨钉固接的标记物，通过捕捉标记物的运动来分析关节运动，此方法虽然降低了皮肤滑移影响，但具有严重侵入性和创伤性，对于正常受试者较难完成相应测试；近年来，计算机断层扫描(ct)和核磁共振成像(mri)和x光等医学成像技术被引入到活体肌肉骨骼关节生物力学研究中，无需使用侵入性标记即可对在体骨骼的形态和位置进行三维量化；但是这些技术多数局限于静态或准静态、二维等限制。
5.近年来，双平面动态x光成像系统(dual fluoroscopic imaging system，dfis)在生物力学研究中获得了发展和应用；该系统是使用两台正交x光机同步捕捉目标关节运动，得到运动关节的双平面动态x光图像，使用3d-2d配准技术，分析骨骼模型在三维空间中的运动；所谓3d-2d配准，是指在虚拟三维空间中重建双平面x光系统环境，包括光源位置及影像增强器位置，基于骨骼轮廓或图像灰度，将三维骨模型与每帧双平面x光图像能正确地匹配，也称为基于模型的追踪。
6.dfis能满足在体、无创、三维、运动测试，分析关节三维运动精度较高。因此，dfis成为了生物运动测试领域的研究热点；但是，此方法的关键技术是3d-2d配准，现有的配准技术多是借助三维软件rhino、maya等构建虚拟双平面环境，进行手动逐帧匹配，需在每一帧的x光图像手动绘制骨骼轮廓，然后移动或旋转骨骼模型，使之轮廓与图像轮廓对齐，而
dfis采集频率较高，在特定动作下捕捉到的x光图像可达几百帧甚至上千帧，这就导致了手动配准耗时费力，且误差较大，是当前使用dfis研究关节运动的技术瓶颈。
7.此外，对于完成3d-2d配准后，一般得到的是骨骼的空间6自由度数据，并不能反映关节(相邻的两个骨骼)的相对运动，需要借助数学分析方法，完成关节三维运动分析；对于关节三维运动的表征与计算，不同的研究采用不同的方法，如jcs(joint coordinate systems)、helical angles等，使得不同的研究结果之间的可比性较差。
8.综观上述目前3d-2d配准技术和关节三维运动分析现状，急需一种可以快速高效的半自动配准方法和具有生物力学意义的关节三维运动求解算法。


技术实现要素：

9.本发明的目的是提供一种基于双平面x光的半自动配准方法，解决目前3d-2d配准多为手动操作而导致的效率低、误差大的问题，并在此基础上，开发了基于双平面x光半自动配准数据的关节三维运动求解算法，以精确计算关节三维六自由度运动，分析关节的三维运动特征。
10.一种基于双平面x光的半自动配准方法和关节三维运动求解算法，所述的基于双平面x光的半自动配准方法包括以下步骤：
11.s1：ct图像处理，包括：
12.s11：创建目标区域：基于目标关节的ct图像，在解剖学平面(矢状面、额状面或横断面)针对某个骨骼创建骨骼外轮廓，逐帧描绘截面轮廓，或隔帧描绘后使用插值法自动填充漏掉的帧，直至在解剖学平面把目标骨骼所包含的所有帧都创建外轮廓，骨骼截面外轮廓所包含的区域即为目标区域；
13.s12：更改像素：将创建目标区域的所有帧，更改目标区域之外的像素值，使之变为透明色，目标区域内像素值保持不变；
14.s13：合成三维体：将更改像素之后的所有帧依据ct图像的体素值合成三维空间体积(长方体)；
15.s14：自定义解剖学坐标系：由s11所述的目标区域重建三维实体骨骼模型，对目标关节所包含的实体骨骼模型建立解剖学坐标系，对于人体典型的关节，如肩关节、髋关节、膝关节和踝关节复合体等，可参照国际生物力学学会(isb)所推荐的标准方法；
16.将目标关节所包含的每块骨骼做s1所述的图像处理操作。
17.s2：双平面x光图像的处理，包括：
18.s21：图像拆分：采集的目标关节动态x光图像一般为图像集合或视频，将图像集合或视频拆分为图像序列；
19.s22：畸变矫正：试验之前，空间标定采集孔阵网格的x光图像，基于孔阵网格文件数据，对原始双平面x光图像序列进行畸变矫正，失真校正算法将孔阵网格图像的孔之间的间距与已知的实际间距值进行比较，并计算用于校正图像的变换矩阵，然后将该矩阵应用于收集到的所有图像序列。
20.s3：重建实验室三维环境，包括：
21.s31：建立全局坐标系：实验前的标定使用内嵌64个金属钢珠和4个特殊形状金属标记物的乐高立方体，将标定立方体中心作为全局实验室坐标系原点，x轴平行人体上下方
向，y轴平行人体行走前后方向，z轴平行人体内外侧方向；
22.s32：建立局部坐标系：由s13所述的将更改像素之后的所有帧合成三维空间体(长方体)，将长方体置于重建的实验室三维环境中，长方体的中心定义为局部坐标系的原点，局部坐标系的x、y和z轴的方向平行于长方体的长、宽和高所在的直线。在配准时，随着x光图像帧数的改变，局部坐标系的位姿也不断地改变；
23.s33：确定光源和双平面图像位置：乐高立方体内嵌64个金属钢珠和4个特殊形状金属标记物，对立方体x射线图像中的钢珠和标记物的中心点进行数字化处理，并与中心点的已知位置进行比较，采用直接线性变换(dlt)定义三维物体空间与二维图像平面之间的线性变换，从而确定相机(光源)和成像平面(双平面图像)的位姿矩阵，即在全局环境下确定光源和双平面图像位置；
24.在重建的实验室三维环境中，由光源发射虚拟x光，双平面图像接收x光，骨骼所在的三维体在虚拟x光的作用下只对骨骼可视化，骨骼之外的体积透明化，以此产生虚拟骨模型。
25.s4：半自动3d-2d配准，包括：
26.s41：添加边缘检测和对比度滤波器：为了提高配准精度，将边缘检测滤波器和对比度增强滤波器应用于骨骼模型和双平面图像序列，可以增强x光图像中骨骼边缘强度和图像清晰度，以提高对比效果。
27.s42：手动配准初始帧，计算相似度：在重建的实验室三维环境中，选取骨骼轮廓和灰度相对清晰的一帧双平面x光图像作为初始帧，通过平移和旋转虚拟骨模型，使得在两个二维视图中(从光源方向，x光图像的正视图)，虚拟骨模型的轮廓与两个x光图像中的骨骼轮廓完全重合，x光图像中骨骼不同的区域灰度不同，在不同视角下，虚拟骨模型的不同区域灰度也不相同，基于骨骼轮廓和灰度，使用下坡单纯形优化算法或粒子群优化算法计算当前帧的虚拟骨模型和x光图像的匹配相似度，设定合格相似度值，当计算出的相似度值小于合格相似度值，即完成初始帧的配准，否则重新调整虚拟骨模型的空间位姿，再进行计算匹配相似度，直至计算出的相似度值小于合格相似度值，记录初始帧对应的局部坐标系的位姿；
28.s43：自动配准整个图像序列：基于配准的初始帧，分阶段对图像序列执行优化算法，例如自动配准接下来的5至10帧，随着帧数改变，自动调整虚拟骨模型的位姿，使其与两个x光图像中的骨骼轮廓完全重合，且灰度强弱匹配，自动记录当前帧对应的局部坐标系的位姿，并在每一帧计算匹配相似度，若匹配相似度小于合格相似度值，即完成此阶段的配准，如果谋帧匹配相似度不小于合格相似度值，则需手动调整虚拟骨模型的位姿，重新计算匹配相似度，直至此阶段的每一帧图像对应的相似度值都小于合格相似度值。分阶段执行自动配准直至完成整个图像序列的配准，输出每帧对应的虚拟骨模型的位姿矩阵；
29.将目标关节所包含的每块骨骼均需完成s4所述的半自动配准；
30.虚拟骨模型的位姿矩阵是基于局部坐标系和全局坐标系之间的位置关系得到的，其并不能反映关节的运动特征，为了使计算出的关节角和关节位移具有运动学意义和生物力学意义，需要一系列的坐标变换，将自定义解剖学坐标系和双平面环境的全局坐标系相关联。
31.所述的关节三维运动求解算法包括以下步骤：
32.s5：坐标变换，包括：
33.s51：解剖学坐标系/ct坐标系：目标关节对应的两个骨骼分别记为a和b，自定义解剖学坐标系是在ct坐标系下定义的，首先确立骨骼a和骨骼b的解剖学坐标系a
an
和b
an
与ct扫描坐标系(ct)之间的转换矩阵。对于骨骼a和骨骼b，从解剖学坐标系到ct坐标系的转换关系表示为：
34.a
an
/ct
35.b
an
/ct
36.s52：ct坐标系/局部坐标系：骨骼a和骨骼b所对应的局部坐标系分别记为la和lb，从ct坐标系到骨骼a和骨骼b所对应的局部坐标系的转换关系表示为：
37.ct/la38.ct/lb39.s53：局部坐标系/全局坐标系：由s4所述的半自动3d-2d配准，输出的是从局部坐标系到全局坐标系的转换矩阵，骨骼a和骨骼b所对应的局部坐标系和全局坐标系(gc)之间的关系表示为：
40.la/gc
41.lb/gc
42.将s51、s52和s53所述的变换矩阵组合以将解剖学坐标系和全局坐标系相关联，骨骼a和骨骼b所对应的转换关系分别为：
43.a
an
/gc＝a
an
/ct
×
ct/la×
la/gc
44.b
an
/gc＝b
an
/ct
×
ct/lb×
lb/gc
45.s6：计算关节角和关节位移，包括：
46.s61：提取欧拉角：构成目标关节的骨骼a和骨骼b一般指近端骨骼和远端骨骼，关节角度是通过从近端到远端的整体转换来确定的，整体转换矩阵记为r
joint
，近端骨骼解剖学坐标系和全局坐标系之间的转换矩阵记为r
proxima1
，远端骨骼解剖学坐标系和全局坐标系之间的转换矩阵记为r
distal
；
[0047][0048]
绕三个坐标轴的旋转角记为：θ
x
、θy和θz，根据转换矩阵r
joint
各元素之间的关系，计算欧拉角θ
x
、θy和θz；
[0049]
s62：计算关节位移：骨骼a的解剖学坐标系原点记为o1，骨骼b的解剖学坐标系原点记为o2，关节位移是通过解剖学坐标系原点转换为全局坐标系来确定的，o1和o2在全局坐标系下的三维坐标记为o
1gc
和o
2gc
；
[0050]o1gc
＝gc/la×
la/ct
×
o1[0051]o2gc
＝gc/lb×
lb/ct
×
o2[0052]
将骨骼a和骨骼b解剖学坐标系原点之间的向量投影到近端骨骼解剖学坐标系的每个轴线上，以量化平移，结果表述为在近端的骨的解剖学坐标系三个坐标轴方向上的移动。
[0053]
本发明具有以下有益效果：
[0054]
本发明提出的基于双平面x光的半自动配准方法，相比现有的手动配准方法，大大
提高了处理效率，节省了人力和时间，有效提高了配准精度，为基于双平面动态x光成像系统的生物三维运动图像处理提供了重要方法和技术参考。
[0055]
本发明所提出的关节三维运动求解算法，使得基于双平面配准的骨骼三维运动数据转化为关节的六自由度相对运动数据，使得研究人员能够完成具有生物力学意义的关节三维运动分析。
附图说明
[0056]
图1为本发明实施例的流程示意图。
[0057]
图2为本发明实施例的ct图像处理过程示意图。
[0058]
图3为本发明实施例的重建实验室三维环境示意图。
[0059]
图4为本发明实施例的半自动3d-2d配准具体操作过程示意图。
具体实施方式
[0060]
如图1、图2、图3和图4所示，一种基于双平面x光的半自动配准方法和关节三维运动求解算法，所述的基于双平面x光的半自动配准方法包括以下步骤：
[0061]
s1：ct图像处理，如图2所示，包括：
[0062]
s11：创建目标区域：图2以踝关节为例，采集踝关节的ct图像，在解剖学平面(矢状面、额状面或横断面)对胫骨创建骨骼外轮廓，逐帧描绘截面轮廓，或隔帧描绘后使用插值法自动填充漏掉的帧，直至在解剖学平面把骨骼所包含的所有帧都创建外轮廓。胫骨截面外轮廓所包含的区域即为目标区域；
[0063]
s12：更改像素。将创建目标区域的所有帧，更改目标区域之外的像素值，使之变为透明色，目标区域内像素值保持不变；
[0064]
s13：合成三维体：将更改像素之后的所有帧依据ct图像的体素值合成三维空间体积(长方体)；
[0065]
s14：自定义解剖学坐标系：由s11所述的目标区域重建三维实体骨骼模型，对目标关节所包含的实体骨骼模型建立解剖学坐标系，对于人体典型的关节，如肩关节、髋关节、膝关节和踝关节复合体等，可参照国际生物力学学会(isb)所推荐的标准方法；
[0066]
将目标关节所包含的每块骨骼做s1所述的图像处理操作。
[0067]
s2：双平面x光图像的处理，包括：
[0068]
s21：图像拆分：采集的目标关节动态x光图像一般为图像集合或视频，将图像集合或视频拆分为图像序列；
[0069]
s22：畸变矫正：试验之前，空间标定采集孔阵网格的x光图像，基于孔阵网格文件数据，对原始双平面x光图像序列进行畸变矫正，失真校正算法将孔阵网格图像的孔之间的间距与已知的实际间距值进行比较，并计算用于校正图像的变换矩阵，然后将该矩阵应用于收集到的所有图像序列。
[0070]
s3：重建实验室三维环境，如图3所示，包括：
[0071]
s31：建立全局坐标系：实验前的标定使用内嵌64个金属钢珠和4个特殊形状金属标记物的乐高立方体，将标定立方体中心作为全局实验室坐标系原点，x轴平行人体上下方向，y轴平行人体行走前后方向，z轴平行人体内外侧方向；如图3所示的在重建实验室三维
环境中的全局坐标系；
[0072]
s32：建立局部坐标系：由s13所述的将更改像素之后的所有帧合成三维空间体(长方体)，将长方体置于重建的实验室三维环境中，长方体的中心定义为局部坐标系的原点，局部坐标系的x、y和z轴的方向平行于长方体的长、宽和高所在的直线；如图3所示的在重建实验室三维环境中的两个骨骼的局部坐标系分别为la和lb。在配准时，随着x光图像帧数的改变，局部坐标系的位姿也不断地改变；
[0073]
s33：确定光源和双平面图像位置：乐高立方体内嵌64个金属钢珠和4个特殊形状金属标记物，对立方体x射线图像中的钢珠和标记物的中心点进行数字化处理，并与中心点的已知位置进行比较，采用直接线性变换(dlt)定义三维物体空间与二维图像平面之间的线性变换，从而确定相机(光源)和成像平面(双平面图像)的位姿矩阵，即在全局环境下确定光源和双平面图像位置；图3所示的在三维环境中的光源和双平面图像的相对位置；
[0074]
在重建的实验室三维环境中，由光源发射虚拟x光，双平面图像接收x光，骨骼所在的三维体在虚拟x光的作用下只对骨骼可视化，骨骼之外的体积透明化，以此产生虚拟骨模型，如图3所示的在重建实验室三维环境中对胫骨和距骨可视化，产生胫距关节虚拟骨模型。
[0075]
s4：半自动3d-2d配准，如图4所示步骤，包括：
[0076]
s41：添加边缘检测和对比度滤波器：为了提高配准精度，将边缘检测滤波器和对比度增强滤波器应用于虚拟骨模型和双平面图像序列，可以增强x光图像中骨骼边缘强度和图像清晰度，以提高对比效果；
[0077]
s42：手动配准初始帧，计算相似度：在重建的实验室三维环境中，选取骨骼轮廓和灰度相对清晰的一帧双平面x光图像作为初始帧，如图4选取踝关节中立位作为初始帧，通过平移和旋转虚拟骨模型，使得在两个二维视图中(从光源方向，x光图像的正视图)，虚拟骨模型的轮廓与两个x光图像中的骨骼轮廓完全重合，x光图像中骨骼不同的区域灰度不同，在不同视角下，虚拟骨模型的不同区域灰度也不相同，基于骨骼轮廓和灰度，使用下坡单纯形优化算法或粒子群优化算法计算当前帧的虚拟骨模型和x光图像的匹配相似度，设定合格相似度值，当计算出的相似度值小于合格相似度值，即完成初始帧的配准，否则重新调整虚拟骨模型的空间位姿，再进行计算匹配相似度，直至计算出的相似度值小于合格相似度值，记录初始帧对应的局部坐标系的位姿；
[0078]
s43：自动配准整个图像序列：基于配准的初始帧，分阶段对图像序列执行优化算法，例如自动配准接下来的5至10帧，随着帧数改变，自动调整虚拟骨模型的位姿，使其与两个x光图像中的骨骼轮廓完全重合，且灰度强弱匹配，自动记录当前帧对应的局部坐标系的位姿，并在每一帧计算匹配相似度，若匹配相似度小于合格相似度值，即完成此阶段的配准，如果谋帧匹配相似度不小于合格相似度值，则需手动调整虚拟骨模型的位姿，重新计算匹配相似度，直至此阶段的每一帧图像对应的相似度值都小于合格相似度值，分阶段执行自动配准直至完成整个图像序列的配准，输出每帧对应的虚拟骨模型的位姿矩阵；
[0079]
将目标关节所包含的每个骨骼均需完成s4所述的半自动配准；
[0080]
虚拟骨模型的位姿矩阵是基于局部坐标系和全局坐标系之间的位置关系得到的，其并不能反映关节的运动特征，为了使计算出的关节角和关节位移具有运动学意义和生物力学意义，需要一系列的坐标变换，将自定义解剖学坐标系和双平面环境的全局坐标系相
关联。
[0081]
所述的关节三维运动求解算法包括以下步骤：
[0082]
s5：坐标变换，包括：
[0083]
s51：解剖学坐标系/ct坐标系：目标关节对应的两个骨骼分别记为a和b，自定义解剖学坐标系是在ct坐标系下定义的，首先确立骨骼a和骨骼b的解剖学坐标系a
an
和b
an
与ct扫描坐标系(ct)之间的转换矩阵；对于骨骼a和骨骼b，从解剖学坐标系到ct坐标系的转换关系表示为：
[0084]aan
/ct
[0085]ban
/ct
[0086]
s52：ct坐标系/局部坐标系：骨骼a和骨骼b所对应的局部坐标系分别记为la和lb，从ct坐标系到骨骼a和骨骼b所对应的局部坐标系的转换关系表示为：
[0087]
ct/la[0088]
ct/lb[0089]
s53：局部坐标系/全局坐标系：由s4所述的半自动3d-2d配准，输出的是从局部坐标系到全局坐标系的转换矩阵，骨骼a和骨骼b所对应的局部坐标系和全局坐标系(gc)之间的关系表示为：
[0090]
la/gc
[0091]
lb/gc
[0092]
将s51、s52和s53所述的变换矩阵组合以将解剖学坐标系和全局坐标系相关联，骨骼a和骨骼b所对应的转换关系分别为：
[0093]aan
/gc＝a
an
/ct
×
ct/la×
la/gc
[0094]ban
/gc＝b
an
/ct
×
ct/lb×
lb/gc
[0095]
s6：计算关节角和关节位移，包括：
[0096]
s61：提取欧拉角：构成目标关节的骨骼a和骨骼b一般指近端骨骼和远端骨骼，关节角度是通过从近端到远端的整体转换来确定的，整体转换矩阵记为r
joint
，近端骨骼解剖学坐标系和全局坐标系之间的转换矩阵记为r
proximal
，远端骨骼解剖学坐标系和全局坐标系之间的转换矩阵记为r
distal
；
[0097][0098]
绕三个坐标轴的旋转角记为：θ
x
、θy和θz，根据转换矩阵r
joint
各元素之间的关系，计算欧拉角θ
x
、θy和θz；
[0099]
s62：计算关节位移：骨骼a的解剖学坐标系原点记为o1，骨骼b的解剖学坐标系原点记为o2，关节位移是通过解剖学坐标系原点转换为全局坐标系来确定的，o1和o2在全局坐标系下的三维坐标记为o
1gc
和o
2gc
：
[0100]o1gc
＝gc/la×
la/ct
×
o1[0101]o2gc
＝gc/lb×
lb/ct
×
o2[0102]
将骨骼a和骨骼b解剖学坐标系原点之间的向量投影到近端骨骼解剖学坐标系的每个轴线上，以量化平移，结果表述为在近端的骨的解剖学坐标系三个坐标轴方向上的移动。
[0103]
本发明提出的基于双平面x光的半自动配准方法，相比现有的手动配准方法，大大提高了处理效率，节省了人力和时间，有效提高了配准精度，为基于双平面动态x光成像系统的生物三维运动图像处理提供了重要方法和技术参考。
[0104]
本发明所提出的关节三维运动求解算法，使得基于双平面配准的骨骼三维运动数据转化为关节的六自由度相对运动数据。使得研究人员能够完成具有生物力学意义的关节三维运动分析。