2. 西北工业大学 软件与微电子学院, 陕西 西安 710072;
3. 西北工业大学 管理学院, 陕西 西安 710072;
4. 西安电子科技大学 计算机学院, 陕西 西安 710071;
5. 西安点云生物科技有限公司, 陕西 西安 710077
由于现代人的生活工作等原因, 脊椎病越来越常见和高发。近年来, 随着计算机可视化在医学影像上的不断应用, 脊椎外科诊疗已经有了很大的提高。但是由于人体脊椎生理结构固有的复杂性, 如切面多、镂空多等, 脊椎外科手术需要复杂且精确的骨定位, 外科手术的难度依然不小。虽然国内外在医学图像三维可视化、医学三维图像的分割等方向进行了大量且卓有成效的研究, 但在关于人体脊椎三维模型特征点的标注方面研究仍然相对较少[1-10]。孟勃等提出的基于改进的尺度不变特征变换特征点匹配方法保持了电子图像的配准精度[1-3], Ballerini等人提出了通过提取模型轮廓曲线的方法, 从而间接得到模型的采样特征点[4-6], 李康等[7]提出的提出了一种基于深度图像的人脸三维模型特征点标定方法, 该方法首先生成人脸三维模型的二维深度图像, 然后采用SUSAN算子在该图像上标定特征点; 翟凤文等[8]提出的置信传播与特征点形状特征相结合去除SIFT特征点匹配错误的方法基于三维SIFT特征匹配方法进行改进, 但上述方法基本是基于模型配准的准则, 如果模型数据较大, 则算法要求计算量也随之变大, 根本满足不了医学研究领域对实时性的要求; 并且由于医学个体的差异性, 完全自动提取方法也存在特征点自动提取后坐标差异较大等根本缺点, 无法满足临床医学的需要; 韩霜等[10]提出基于高斯曲率的半自动特征点标注方法, 方法虽然能够对标注结果自适应调节, 但高斯曲率却不能极大表述模型几何结构, 故不能准确反映模型的局部几何特征, 且精度不高。
为了帮助医生更加准确地诊断、治疗, 给临床提供更为准确地科学度量数据, 本文提出了一种基于法曲率极大值和向量内积的脊椎特征点自动标注方法, 该方法首先通过对脊椎CT数据进行三维重建和模型分割, 从而获得需要进行特征点识别的特定椎体模型; 计算手动选定特征点周围极小半径r球体内各个顶点的高斯曲率和平均曲率流, 从而得到r球体范围内各个顶点的法曲率极大值; 按照极大值曲率大小进行顶点降序排列, 选出法曲率极大值最大的n个候选顶点; 以空间原点为向量的起点, 对候选顶点与手动选取点之间作向量内积, 从而获得两顶点向量之间的几何夹角角度, 以夹角角度最小的候选顶点来替代手动选取的点来作为模型该区域的特征点。
1 基于法曲率极大值和向量内积的人体三维脊椎模型特征点自动识别方法 1.1 特征点标定在三维图像学中, 特征点的定义是指在三维模型中, 用来表征局部邻域关系及特性的, 具有一定几何意义的关键三维坐标点, 特征点矩阵可用来存储或者表示三维模型的几何轮廓、生理特征等属性信息。
本文根据Cootes提出的标记点分类定义, 并借鉴在国际上得到广泛认同的颅脑特征点标注标准MPEG-4(moving picture experts group)的人脸定义参数(FDP)和人体生理学脊椎特征, 标记出脊椎三维特征点。图 1和图 2分别是在骨骼解剖学中的脊椎腰骶段的右侧面观和上面观特征点示意图, 如图 1和图 2所示。
1.2 基于基于法曲率极大值和向量内积的特征点自动识别方法人体脊椎模型特征点特指包含了更多模型局部几何信息的骨点, 这些点具有特征突出、稳定性好等特点, 同时, 也是自然曲率较大的点。这些点由于曲率较大, 三维模型表面在该点处的弯曲程度自然也较大, 即就是说, 该点能最大程度的表现三维模型的大致轮廓和局部模型变化, 从而准确反映该点的特征情况。但采用纯手工的方法标注效率不高, 准确性也略差; 完全自动提取方法也存在坐标提取差异大等缺点。为了解决这些问题, 本节在手动拾取特征点的基础上提出了一种基于曲率极大值和向量内积的自动调整方法, 以保证拾取特征点的准确性。方法通过计算手动选定模型特征点周围极小半径r球体内各个顶点的高斯曲率和平均曲率流, 从而得到r球体范围内各个顶点的法曲率极大值; 然后按照极大值曲率大小进行顶点降序排列, 选出法曲率极大值最大的n个候选顶点; 对候选顶点与手动选取点一一作向量内积, 从而获得两顶点向量之间的几何夹角角度, 以夹角角度最小的候选顶点来替代手动选取的点来作为模型该区域的特征点, 由于在极小r半径内, 向量之间的夹角越小, 则表示两顶点之间逼近程度越高, 从而可以最大程度保证该点特征变化情况的准确性。
1.2.1 三角网格表面的极大值曲率计算在计算机三维图像处理中, 物体的曲面表结构都是由一簇网格的线性逼近值所组成的。因此, 可以考虑将三角面片网格的每个顶点的度量性质看做是此点局部邻域的平均度量。
如考虑计算图 3中点xi邻域的极大曲率值, 用冯洛诺伊图(Voronoi diagram)内有且仅含有一个的离散点极大值曲率来做替换。
假设kmax和kmin分别为曲面上一点处法曲率的最大值、最小值, 则根据曲率的定义, 有:
(1) |
(2) |
联立(1)式和(2)式, 则可以求得(3)式和(4)式:
(3) |
(4) |
式中,kG为高斯曲率值, kH为平均曲率值。
1.2.2 向量内积几何夹角的计算利用点积公式(5)计算手动选取点与候选点的2个向量内积:
(5) |
式中, O为坐标系原点, p0为手动选取点, pi为某个法曲率极大值最大的候选点;
利用公式(6)计算两向量的夹角θ:
(6) |
从而得到, 候选点与手动选取点之间的几何夹角角度, 其值大小决定了特征点的准确度, 值越小, 表明方法拾取的点更加集中, 即更加趋近于真实特征点, 准确度也越高。
当有极小概率事件-几何夹角角度相同时, 则根据空间距离公式(7)来进行判定, 欧几里得距离较小的则为最优点。
(7) |
本标注改进方法的主要步骤如下:手动选择模型某处特征点范围内的任一点p0, 以p0为中心, 给定极小r为半径, 计算在该球体范围内的所有点的高斯曲率值KG和平均曲率值KH, 从而可以计算得到每个点的法曲率极大值Kmax, 一般情况下, 曲率较大, 曲线拐点信息越明确, 几何特征愈发明显, 对法曲率极大值Kmax进行降序排序。选取Kmax最大的n(n<=3)个点, 对这n个点, 分别与手动选取点p0作内积, 夹角最小的候选点即为所求的点, 可以替换p0作为新的模型特征点。
Step1 获得脊椎三维图像表面的结构信息facets。
Step2 定义的各个变量曲率KG、KH和领域面积AM, 边向量长度Length, 角度向量Angle并初始化。
Step3 计算领域面积AM。通过遍历, 得到模型的每一个三角面片的坐标点向量, 从而可以计算出边向量以及三角面片内角度数和, 并调用数学计算库函数三角面片求积得到AM。
Step4 遍历计算模型上每个点的高斯曲率值KG。
Step5 计算向量(e1, e2)、(e2, e0)和(e0, e1)的边向量长度Length
Step6 计算三角面片P(x, y, z)的法向量n_f, 及其邻域三角面片的法向量n_n, 依据点积公式cos=n_f * n_n从而可以计算出n_f*n_n的点积cos, 并且依据公式sin =(n_f * n_n) * P′(x, y, z), 可得到sin值。调用atan2(sin, cos)函数, 获得角度向量Angle。
Step7 遍历计算模型上每个点的平均曲率值KH。
Step8 利用公式(3)遍历计算得到模型上每个点的法曲率极大值Kmax。
Step9 通过每个顶点的ID值可以将curvature数组中存放的极大法曲率值与脊椎三维表面上的特征点对应起来。计算r范围内所有ID的法曲率极大值Kmax, 通过降序排序, 找出该球体范围内法曲率极大值最大的n(n<=3)个点。
Step10 利用公式(5)将选出的法曲率极大值最大的n(n<=3)个点与手动选取点做向量
内积, 并由公式(6)分别计算与手动选取点的向量夹角, 夹角最小的候选点, 即为经过方法修正的特征点。
2 方法验证与分析本方法验证采用主频为2.60GHZ的4GB四核AMD主机, 操作系统为WINDOSW 7, 采用Microsoft MFC library作为GUI, 基于VC++开发语言, 并利用开源的VTK开发包作为图像显示增强库。利用以上技术, 我们开发了专门的交互软件用以验证该改进方法。
2.1 改进后的特征点拾取结果在VC++平台下, 图 4是利用上述改进后的方法程序对脊椎三维面模型进行特征点拾取的效果图:
图 4中圆圈内的拾取点为1号特征点。左下角显示的拾取数据里, 前2行为纯手工选取的特征点的坐标及其极大值曲率。后2行为经过本方法调整过之后选取的特征点的坐标、极大值曲率和角度值。
2.2 改进前后结果的对比分析为了验证本方法的准确性, 做出如下对比实验。特意选取如下3种特征点:1号特征点(模型基准点)、7号特征点(特征较不明显的点)和32号特征(模型轮廓特征点)点, 在不同绘制方位的条件下, 采用改进前的手动拾取特征点方法进行特征点标定, 实验效果如图 5所示:
图 5中, 图a)、b)、c)分别是在不同绘制条件下进行了3次实验后的截图。图中显示的数据分别为对1号、7号和32号特征点进行标定的特征点坐标及其曲率。
为了更好地进行对比, 利用改进前手动拾取方法和本方法, 分别对No.1、7、32进行5组标记, 得到标记的三维坐标点和对应的极大值曲率以及利用本方法调整后的特征坐标点和极大值曲率, 具体数值如表 1所示。
特征点序号 | P(x, y, z) | 高斯曲率 | 平均曲率 | 极大值曲率流 | 调整后的P′(x, y, z) | 调整后的极大值曲率流 |
No.1 | P1(68.617 259, 94.271 348, 43.866 595) | 0.092 735 | 0.412 750 | 0.691 369 | P′(69.372 986, 93.936 020, 44.156 120) | 1.096 345 |
P2(67.917 100, 94.060 941, 43.093 403) | 0.002 077 | 0.078 837 | 0.145 094 | P′(69.944 412, 94.316 490, 42.861 012) | 1.388 327 | |
P3(68.266 710, 93.840 342, 44.388 124) | 0.011 477 | 0.048 556 | 0.135 387 | P′(68.947 128, 93.506 936, 44.622 269) | 0.983 049 | |
P4(68.267 095, 94.271 123, 43.868 031) | 0.009 367 | 0.072 781 | 0.167 443 | P′(69.372 986, 93.936 020, 44.156 120) | 1.096 345 | |
P5(68.267 407, 94.512 628, 43.445 553) | 0.007 631 8 | 0.025 400 8 | 0.065 577 2 | P′(68.195 251, 94.839 355, 42.868 443) | 1.693 831 | |
No.7 | P1(58.901 743, 73.082 046, 48.452 738) | 0.112 116 | 0.149 657 | 0.350 621 | P′(59.408 794, 73.907 570, 49.243 549) | 1.331 628 |
P2(59.224 789, 73.062 200, 48.093 790) | 0.068 149 | 0.119 758 | 0.242 043 | P′(59.517 223, 73.743 927, 48.833 157) | 0.724 517 | |
P3(59.385 671, 73.090 835, 47.930 803) | 0.045 982 | 0.233 578 | 0.326 188 | P′(59.517 223, 73.743 927, 48.833 157) | 0.724 517 | |
P4(59.069 080, 72.872 247, 47.837 257) | 0.001 001 | 0.173 883 | 0.350 621 | P′(59.517 223, 73.743 927, 48.833 157) | 0.724 517 | |
P5(58.433 019, 72.536 799, 47.518 840) | 0.057 150 | 0.319 407 | 0.531 236 | P′(59.225 231, 73.472 542, 48.676 968) | 0.717 429 | |
No.32 | P1(44.980 407, 43.096 873, 69.899 599) | 0.072 504 | 0.339 573 | 0.383 041 | P′(44.794 018, 41.827 137, 69.070 534) | 1.492 776 |
P2(45.110 727, 43.105 876, 70.027 456) | 0.029 067 | 0.229 463 | 0.383 041 | P′(44.867 661, 42.867 104, 69.194 344) | 1.459 913 | |
P3(45.406 852, 42.500 475, 70.079 444) | 0.002 911 | 0.075 997 | 0.129 519 | P′(44.794 018, 41.827 137, 69.070 534) | 1.492 776 | |
P4(45.663 887, 42.973 022, 70.350 514) | 0.030 788 | 0.219 997 | 0.352 702 | P′(44.964 773, 42.077 076, 69.538 544) | 1.848 617 | |
P5(45.658 297, 43.277 567, 70.381 229) | 0.175 850 | 0.804 599 | 0.313 514 | P′(44.156 033, 42.608 818, 70.195 366) | 1.884 813 |
由表 1可以看出, 对位于模型表面上的基准点(No.1), 其曲率往往会偏大一些。对位于平滑曲面上的特征较不凸出的点(No.7), 其曲率往往会偏小一些; 针对同一特征点, 其极大值曲率均要大于其高斯曲率、平均曲率, 而调整后的极大值曲率流则要明显大于相对应的手动标注极大值曲率流, 如图 6所示。
同理, 运用文献[10]的方法与本文方法对相同标记点在同等条件下进行5组标记, 其结果如图 7、图 8所示。
在图 7中, 图中显示的数据分别为使用文献[10]的方法对模型1号、7号和32号特征点进行修正标注后的特征点坐标及其曲率。
在图 8中, 图中显示的数据分别为使用本改进方法对模型1号、7号和32号特征点进行调整标注的特征点坐标和调整后的极大值曲率、以及夹角角度。同样, 为了突出本方法的优势, 将文献[10]标注方法与本发明方法在同一实验条件下, 对特征点进行标注, 具体参数结果如表 2和图 9所示。
特征点序号 | P(x, y, z) | 文献[10]调整后的P′(x, y, z) | 文献[10]的方法计算得到的极大值曲率流 | 文献[10]的夹角角度/(°) | 本方法调整后的P′(x, y, z) | 本方法调整后的极大值曲率流 | 本方法调整后的夹角角度/(°) |
No.1 | P1(68.617 259, 94.271 348, 43.866 595) | P′(67.601 565, 94.892 846, 43.853 074) | 1.102 646 | 0.009 319 | P′(69.372 986, 93.936 020, 44.156 120) | 1.096 345 | 0.006 689 |
P2(67.917 100, 94.060 941, 43.093 403) | P′(68.803 757, 94.842 346, 42.855 259) | 1.397 646 | 0.012 733 | P′(69.944 412, 93.316 490, 42.861 012) | 1.388 327 | 0.011 256 | |
P3(68.266 710, 93.840 342, 44.388 124) | P′(69.372 986, 93.936 020, 44.156 120) | 1.001 878 | 0.007 725 | P′(68.947 128, 93.506 936, 44.622 269) | 0.983 049 | 0.005 876 | |
P4(68.267 095, 94.271 123, 43.868 031) | P′(66.944 412, 94.316 490, 42.861 012) | 1.098 327 | 0.010 477 | P′(69.372 986, 93.936 020, 44.156 120) | 1.096 345 | 0.008 817 | |
P5(68.267 407, 94.512 628, 43.445 553) | P′(66.944 677, 94.316 771, 42.861 287) | 1.702 765 | 0.007 973 | P′(68.195 251, 94.839 355, 42.868 443) | 1.693 831 | 0.005 361 | |
No.7 | P1(58.901 743, 73.082 046, 48.452 738) | P′(57.760 895, 73.172 516, 49.520 557) | 1.348 391 | 0.014 807 | P′(59.408 794, 73.907 570, 49.243 549) | 1.331 628 | 0.002 709 |
P2(59.224 789, 73.062 200, 48.093 790, 48.590 868) | P′(58.878 677, 73.417 770, 48.962 761) | 0.734 965 | 0.008 437 | P′(59.517 223, 73.743 927, 48.833 157) | 0.724 517 | 0.003 657 | |
P3(59.385 671, 73.090 835, 47.930 803) | P′(58.878 677, 73.417 770, 48.962 761) | 0.734 965 | 0.010 587 | P′(59.517 223, 73.743 927, 48.833 157) | 0.724 517 | 0.005 812 | |
P4(59.069 080, 72.872 247, 47.837 257) | P′(58.878 677, 73.417 770, 48.962 761) | 0.734 965 | 0.009 375 | P′(59.517 223, 73.743 927, 48.833 157) | 0.724 517 | 0.004 634 | |
P5(58.433 019, 72.536 799, 47.518 840) | P′(59.225 231, 73.472 542, 48.676 968) | 0.717 429 | 0.004 479 | P′(59.234 511, 72.485 362, 49.124 147) | 0.705 282 | 0.003 021 | |
No.32 | P1(44.980 407, 43.096 873, 69.899 599) | P′(44.469 620, 43.914 757, 68.962 631) | 1.492 776 | 0.013 092 | P′(44.794 018, 41.827 137, 69.070 534) | 1.492 402 | 0.008 816 |
P2(45.110 727, 43.105 876, 70.027 456) | P′(44.409 828, 43.943 066, 69.090 309) | 1.476 093 | 0.015 345 | P′(44.867 661, 42.867 104, 69.194 344) | 1.459 913 | 0.006 479 | |
P3(45.406852, 42.500 475, 70.079 444) | P′(45.167 667, 41.388 000, 69.150 368) | 1.509 841 | 0.007 079 | P′(44.794 018, 41.827 137, 69.070 534) | 1.492 776 | 0.000 795 | |
P4(45.663 887, 42.973 022, 70.350 514) | P′(44.627 914, 43.753 708, 69.364 044) | 1.850 723 | 0.014 675 | P′(44.964 773, 42.077 076, 69.538 544) | 1.848 617 | 0.003 713 | |
P5(45.658 297, 43.277 567, 70.381 229) | P′(44.467 278, 43.749 367, 69.313 614) | 1.893 074 | 0.012 892 | P′(44.156 033, 42.608 818, 70.195 366) | 1.884 813 | 0.005 447 |
由表 2和图 9可以看出, 同一模型点的每次标注, 本方法和文献[10]的极大值曲率流相差不大, 但差别较大的方面是, 本方法的几何夹角均要明显小于文献[10]计算的角度夹角。
3 结论从上述数据中可以看到, 对于同一模型特征点, 其法曲率极大值均大于高斯曲率和平均曲率。另外, 比较表 1的法曲率极大值可知, 与手动拾取相比, 本方法所调整的新坐标点, 其法曲率极大值的绝对值均更大。依据论文前述的极大值曲率值的几何意义, 可知, 本方法标注的特征点包含更多的几何信息, 所拾取的特征点更为准确。另外, 通过比较表 2的极大值曲率流和向量夹角, 也可以看到, 对于不同方位的模型特征点, 在基本确保曲率精度的同时, 本方法的几何夹角均要明显小于文献[10]的方法, 依据之前论述的几何夹角意义可知, 向量之间几何夹角越小, 代表 2个点之间欧几里得距离越靠近, 从而表明新方法标注的点更趋于聚集, 更能逼近实际的特征点信息; 尤其是对于特征不明显的点, 比如No.7特征点, 本方法的改进效果最明显, 即本方法能更好的避免因为人为所产生的误差。经过加权计算得知, 本方法的标注准确度提高了约35%, 算法复杂度为O(nlog2n), 如果r范围内的像素点是几何级倍数的话, 使用本方法, 既能保证精准度的同时又可以保证实时性, 从而, 为医学脊椎疾病的诊疗和判断提供较为准确的科学度量数据。
[1] |
孟勃, 韩广良. 基于改进的尺度不变特征变换特征点匹配的电子稳像方法[J]. 计算机应用, 2012, 32(10): 2817-2820.
Meng Bo, Han Guangliang. Electronic Image Stabilization Method Using Improved Scale Invariant Feature Transform[J]. Journal of Computer Applications, 2012, 32(10): 2817-2820. (in Chinese) |
[2] |
李根, 李文辉. 基于思维进化方法的人脸特征点跟踪[J]. 吉林大学学报, 2015, 45: 606-612.
Li Gen, Li Wenhui. Facial Feature Tracking Based on Mind Evolutionary Algorithm[J]. Journal of Jilin University, 2015, 45: 606-612. (in Chinese) |
[3] |
刘志文, 刘定生, 刘鹏. 应用尺度不变特征变换的多源遥感影像特征点匹配[J]. 光学精密工程, 2013, 21(8): 2146-2153.
Liu Zhwen, Liu Dingsheng, Liu Peng. Sift Feature Matching Algorithm of Multi-Source Remote Image[J]. Opt Precision Eng, 2013, 21(8): 2146-2153. (in Chinese) |
[4] | Ballerinil Calistim, Cordon O. Auto-Matic Feature Extraction from 3D Range Images of Skulls[J]. Compuational Forensics, 2011, 5158: 58-69. |
[5] |
王欣, 张明明, 于晓, 等. 应用改进迭代最近点方法的点云数据配准[J]. 光学精密工程, 2012, 20(9): 2068-2077.
Wang Xing, Zhang Mingming, Yu Xiao, et al. Point Cloud Registration Based on Improved Iterative Closest Point Method[J]. Opt Precision Eng, 2012, 20(9): 2068-2077. (in Chinese) |
[6] | Damas S, Cordon O, Santamaria J, et al. Forensic Identification by Computer-Aided Craninfacial Superimposition:A Survery[M]. Acm Computing Surveys(Csur), 2011: 27. |
[7] |
李康, 尚鹏. 基于深度图像的人脸模型特征点自动标定[J]. 计算机科学, 2014, 41: 287-291.
Li Kang, Shang Peng. Automatic Location of Feature Points on Three-Dimensional Facial Model Based on Depth Image[J]. Computer Science, 2014, 41: 287-291. (in Chinese) |
[8] |
翟凤文, 党建武. 形状特征和置信传播在去除Sift特征点错误匹配中的应用[J]. 计算机辅助设计与图形学学报, 2016, 28: 443-449.
Qu Fengwen, Dang Jianwu. Application of Shape Context and Belief Propagation in Removing Sift Mismatches[J]. Journal of Computer-Aided Design & Computer Graphics, 2016, 28: 443-449. DOI:10.3969/j.issn.1003-9775.2016.03.009 (in Chinese) |
[9] |
王澜, 孙博. 一种基于共线特征点的线阵相机内参标定方法[J]. 红外与激光工程, 2015(6): 134-140.
Wang Lan, Sun Bo. A Novel Method for Calibrating Intrinsic Parameters of Linear Array Cameras Based on Collinear Feature Points[J]. Infrared and Laser Engineering, 2015(6): 134-140. (in Chinese) |
[10] |
韩霜. 三维脊椎模型特征点标注的研究与实现[D]. 西安: 西安电子科技大学, 2013 Han Shuang. Research and Implementation of Tagging Three-Dimensional Model of the Spine with Feature Points[D]. Xi'an, Xidian University, 2013(in Chinese) http://cdmd.cnki.com.cn/article/cdmd-10701-1014327772.htm |
[11] | 李瑞祥. 实用人体解剖彩色图谱[M]. 北京: 人民卫生出版社, 2001. |
[12] |
陈维桓. 微分几何[M]. 北京: 北京大学出版社, 2006.
Chen Weiheng. Differential Geometry[M]. Beijing: Peking University Press, 2006. (in Chinese) |
2. School of Software and Microelectronics, Northwestern Polytechnical University, Xi'an 710072, China;
3. School of Management, Northwestern Polytechnical University, Xi'an 710072, China;
4. College of Computer Science, Xidian University Xi'an 710071, China;
5. Xi'an Particle Cloud Biotechnology Co., Ltd, Xi'an 710072, China