一种弯骨骨折自动配准及内固定钢板预弯建模方法与流程

本发明涉及弯骨的截面分割及配准技术领域，尤其涉及一种弯骨骨折自动配准及内固定钢板预弯建模方法。

背景技术

目前骨折手术普遍采用人工复位和伤肢内固定相结合的方法，这种方法存在的问题是创伤大、出血多、容易引发神经血管损伤等并发症。因此，我们可以利用计算机算法对断骨模型进行虚拟拼接，从而在术前得到钢板的各种几何参数。然而目前的断骨虚拟拼接方法仅能应用于长直骨，在弯骨的截面分割及配准领域还没有较为有效的全自动配准方法。并且在基于脊线预配准的肋骨配准方法中，需要手工拾取种子点来进行截面点集的分割，没有实现全自动的截面点集分割和弯骨配准。

技术实现要素：

根据现有技术存在的问题，本发明公开了一种弯骨骨折自动配准及内固定钢板预弯建模方法，具体包括以下步骤：

s1：采用pca算法提取断骨模型点集的主方向并沿主方向选取聚类种子点,采用k-means聚类方法将模型的聚类种子点集进行聚类计算每个聚类的中心点，利用聚类的中心点来拟合bezier曲线获得断骨模型的走势曲线信息；

s2：在断骨模型的两端分别构造胶囊状的映射模型，将断骨模型内的点沿其法向量方向映射到映射模型上；

s3：对s2中构造的映射模型采用optics算法对映射模型的特征点集进行聚类处理，对映射模型的特征点集进行特征匹配实现两个映射模型的粗配准；

s4：对断骨模型截面部分的点集进行提取和分割：将粗配准中的特征点进行区域自生长得到扩展点集，在扩展点集中进行匹配搜索得到截面点集；

s5：根据分割得到的两个截面点集对断骨模型进行精确配准：利用pca算法提取截面点集的主方向并根据主方向进行预配准，利用icp算法进行精配准；

s6：在断骨模型的断裂部位选取控制点，将控制点范围内的三角面片进行加厚处理得到钢板模型。

进一步的，s2中具体采用如下方式：

根据拟合得到的断骨模型的走势曲线计算该曲线两端的切线，分别以四个切线为轴线构造胶囊形状的四个映射模型，由于四个映射模型具有相同的尺寸，分别将映射模型范围内的点集沿其法向量方向映射到映射模型上，得到四个包含映射点集的映射模型。

进一步的，所述s3中两个映射模型的粗配准过程具体采用如下方式：

将两个映射模型之一进行镜面变换、使两个映射模型在坐标系中进行匹配，对映射模型中圆柱的表面部分进行网格划分、判断映射点集中圆柱部分点集的网格归属，将圆柱部分的点集全部划分到网格，统计每个网格中映射点的数目、将数目大于阈值的网格中所有的点提取出来作为映射点集的特征点，采用optics算法对特征点集进行聚类处理提取聚类中心、对聚类中心点进行遍历，计算聚类之间的匹配度、提取匹配度最大的两个聚类点集并获取旋转角度、根据旋转角度将断骨模型进行旋转来达到粗配准。

进一步的，所述s4中断骨模型截面部分的点集提取和分割过程具体采用如下方式：将粗配准中提取到的两个聚类点集反向映射，在断骨模型表面找到与其对应的点集作为截面特征点，将截面特征点集进行区域自生长得到两个扩展点集，根据两点之间的坐标关系和法向量关系进行匹配点的搜索，得到两个截面的点集。

由于采用了上述技术方案，本发明提供的一种弯骨骨折自动配准及内固定钢板预弯建模方法，本方法可以实现弯骨骨折的自动配准，并根据配准之后的弯骨模型拟合钢板，本发明提出的算法可以全自动地进行弯骨配准，并且鲁棒性高，可以适应各种复杂的断面，具有很高的精度，能够在术前拟合得到手术所需要的钢板模型。

附图说明

为了更清楚地说明本申请实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本申请中记载的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明公开的方法的流程图；

图2为本发明中断骨模型曲线拟合效果图；

图3为本发明中对断骨模型进行直线提取示意图；

图4为本发明中k-means聚类效果示意图；

图5为本发明中构造映射模型示意图；

图6为本发明中粗配准效果示意图；

图7为本发明中正角θ的示意图；

图8为本发明中特征点提取的效果示意图；

图9为本发明中optics聚类效果示意图；

图10为本发明中提取聚类的中心点示意图；

图11为本发明中分割得到的截面点集示意图；

图12为本发明中特征点集区域自生长后的效果示意图；

图13为本发明中预配准效果示意图；

图14为本发明中截面点集精配准的结果示意图；

图15为本发明中断骨模型精配准的效果示意图；

图16为本发明中在模型表面选取控制点示意图；

图17为本发明中拟合得到钢板模型示意图。

具体实施方式

为使本发明的技术方案和优点更加清楚，下面结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚完整的描述：

如图1所示的一种弯骨骨折自动配准及内固定钢板预弯建模方法，具体包括以下步骤：

如图2所示：s1：利用pca算法提取断骨模型点集的主方向并沿主方向选取聚类种子点,采用k-means聚类方法将模型的聚类种子点集进行聚类计算每个聚类的中心点，利用聚类的中心点来拟合bezier曲线获得弯骨模型的走势曲线信息。具体方式为：利用pca算法提取断骨模型点集的主方向。提取效果如图3所示，pca是提取点集主方向的成型算法，其具体步骤为：

计算点集中心点坐标：分别计算断骨模型中所有点的x、y、z坐标的平均值，得到断骨点集中心点的坐标p0(x0,y0,z0)。

特征中心化：将断骨点集的空间坐标构建为3×n矩阵n为点的总数量，分别计算矩阵中每个点的坐标与中心点坐标p0的差值，得到更新后的矩阵a，即

计算协方差矩阵：将矩阵a与其转置矩阵a′相乘，得到协方差矩阵m，即m＝aa′

求协方差矩阵m的特征值和特征向量：最大的特征值对应的向量即为该模型点集的主方向，即提取得到的直线方向。

沿直线选取等距的5个种子点，然后以这5个点作为中心点进行k-means聚类。聚类效果如图4所示。k-means是经典的聚类算法，其具体步骤为：

分别以5个种子点作为中心点创建5个聚类，对点集中的所有点，计算得到离它距离最近的种子点，并将该点加入到该种子点所在的聚类。

所有点全部归类完毕后，重新计算5个类的中心点，并以此作为新的种子点进行下一轮迭代。直到达到设定的迭代次数时，停止迭代。

利用聚类的中心点来拟合bezier曲线。计算曲线表达式，n阶bezier曲线的参数公式为：其中t为参数，取值范围为[0,1)，n为bezier曲线的阶数，pi代表第i个中心点的坐标。

进行插值，在曲线上进行插值，即通过给定不同的t，来计算得到曲线上多个点的坐标，从而近似拟合得到完整的曲线。最终效果如图2所示。

s2：在断骨模型的两端分别构造胶囊状的映射模型，将断骨模型内的点沿其法向量方向映射到映射模型上。具体方式为：计算曲线两端的切线，为方便对每个映射模型进行构建，在每个模型构建时，都要将断骨模型进行移动，使该曲线端点处于坐标原点，切线方向为x轴正方向。

构造的映射模型形状类似于胶囊，中间为圆柱状，两端为半球状，模型的几何中心为断骨模型曲线的端点，模型的轴线为上述计算得到的切线。设置映射模型圆柱部分的横切圆半径为r，圆柱部分的长度为l，则映射模型的坐标方程为

将点集进行映射。对模型范围内的点集，分别沿其所在位置处法向量的方向进行映射，即计算法向量所在直线与映射模型的交点，作为该点的映射点，重复上述步骤，共构建得到4个映射模型，如图5。

s3：对s2中构造的映射模型采用optics算法对映射模型的特征点集进行聚类处理，对映射模型的特征点集进行特征匹配实现两个映射模型的粗配准。具体过程是：将两个映射模型之一进行镜面变换、使两个映射模型在坐标系中进行匹配，对映射模型中圆柱的表面部分进行网格划分、判断映射点集中圆柱部分点集的网格归属、将圆柱部分的点集全部划分到网格，统计每个网格中映射点的数目，将数目大于阈值的网格中所有的点提取出来作为映射点集的特征点，应用optics算法对特征点集进行聚类处理，提取聚类中心、对聚类中心点进行遍历，计算聚类之间的匹配度、提取匹配度最大的两个聚类并获取旋转角度、根据旋转角度将断骨模型进行旋转来达到粗配准。具体算法包括4个步骤：

s31：镜像处理。因为两个相匹配的模型在初始构造时是镜像对称的，所以为了方便后续匹配，需要先做镜像处理。在本方法中，由于映射模型的几何中心都在原点，所以镜像处理时只需将映射点集中每个点的x轴的坐标变为相反数。

s32：特征点提取。对映射模型中圆柱的表面部分进行网格划分，沿圆柱的高(平行x轴方向)纵向划分为nx个区间，沿侧表面圆周划分为ny个区间，所以总的网格数为ngrid＝nx×ny。然后判断映射点集中圆柱部分的网格归属，设映射点的空间坐标为(x,y,z)，则其所属网格的编号为l为圆柱的长度，θ为以y轴正方向为始边，以该点与横切面圆心连线为终边的正角，如图7。将圆柱部分的点集全部划分到网格之后，统计每个网格中映射点的数目，将数目大于阈值(设为3)的网格中所有的点提取出来，作为映射点集的特征点，提取效果如图8所示。

s33：optics聚类。应用optics算法对特征点集进行聚类。具体步骤

(1)遍历特征点集，寻找每个点的临近点(与该点距离小于阈值ε的点)，统计每个点对应的临近点的数目，如果大于阈值m,则该点为核心点，并计算该点的核心距离，p点的核心距离定义为

其中，nε(p)表示点p的临近点集合，表示nε(p)中与p点第i近邻的点。

(2)建立两个队列，分别是有序队列和结果队列。遍历点集集合中的所有点，选择一个未被处理的核心点加入到结果队列中，并将该点的所有临近点加入到有序队列中。计算有序队列中每个点从核心点出发的可达距离，从p0到p1的可达距离的定义为

将有序队列中的点按照可达距离升序排列，如果有序队列不为空，则对其中的点进行遍历，每次取出一个点，若该点是核心点且不在结果队列中，则将该点加入到结果队列，并将其临近点加入到有序队列中，如果其临近点已经在有序队列中且可达距离较小，则更新可达距离，最后对有序队列重新排序。重复上述过程，直至所有的特征点遍历完成。

(3)输出聚类结果。首先设置阈值半径ε′，从结果队列中顺序取出点，如果该点的可达距离不大于ε′，则该点属于当前类别；如果该点的可达距离大于ε′且该点为核心点，则开辟一个新的聚类；如果该点的可达距离大于ε′且不为核心点，则该点为噪音点。聚类的结果如图9所示。

s34：提取聚类中心。首先采用层次聚类的方法将optics聚类结果进行合并处理，合并距离较近的类别，并且删除元素数目少的类别。具体方法是：按编号顺序遍历optics聚类的中心点，对每个类别的中心点，遍历其后续类别的中心点，若两个中心点的距离小于阈值(设为4.0)，则将这两个类进行合并。聚类完成后，遍历并去掉元素数目小于阈值(设为20)的类别，最后对更新完成的聚类重新编号。层次聚类完成之后，计算每个类别的坐标中心点，结果如图10所示。

s35：对映射点集进行特征匹配。具体步骤为：

(1)聚类中心点所在的横切圆的圆心到该中心点的向量为将绕x轴逆时针旋转至z轴正方向经过的角度为α，则对每个聚类的中心点计算其对应的α，将所有的α值进行排序并去掉耦合，加入到角度队列q中。

(2)对其中两个映射模型的特征点集进行比较时，先根据映射模型头部(端点球面部分)的点集数目来判断为何种断面类型。若两个映射模型的头部点集数目均大于阈值(设为300)时，执行如下方法：遍历映射模型二中的头部点集，对每个点pi，在映射模型一的头部点集中寻找其最邻近的点并计算最近距离si，si为球面上的最短距离，并非空间最短距离。引入匹配距离s来表征两个映射点集的匹配程度，s越小代表匹配程度越高。在研究头部点集的方法中，其中nh为映射点集二中点的数目。若两个映射模型的头部点集不都大于阈值时，则分析映射模型的圆柱部分，方法如下：两个映射模型的聚类中心点角度队列分别为q1和q2，对q1q2中的所有角度组合进行遍历，假设当前遍历到的角度组合为(αi,βj),其中αi∈q1，βj∈q2，则将映射模型一逆时针旋转αi，将映射模型二逆时针旋转βj，并将αi，βj对应的中心点所在聚类c,c′设为当前聚类。若c和c′中的点集数目均超过阈值(设为500)，则遍历映射模型二中的圆柱部分点集，对每个点pi′，在映射模型一的圆柱部分点集中寻找其邻域点，可以利用之前步骤中建立的网格模型来提高检索效率。对每种角度组合，分别计算其对应的匹配距离s。首先要判断一种c和c′高度匹配的情况，计算公式为其中a1代表c中匹配点在c′中的点的数目，b1代表c中有匹配点的所有点集数目，a2b2同理。若s′<0.1，则s＝s′。否则的话

m1和m2分别表示c和c′中点的总数目。

(3)对所有的角度组合进行遍历并计算匹配距离s，最小的s对应的角度组合即为该两个映射模型的最优配准角度。对两个断骨模型的共4个映射点集进行两两匹配分析，共进行4组匹配分析，其中匹配距离最小的两个映射点集即为两个断面的映射点集，其对应的配准角度即为粗配准的角度。

(4)将断骨模型按照步骤(3)中求得的角度绕x轴进行逆时针旋转，然后适当平移使粗配准结果更加精确。平移的距离由两个映射模型的头部点集数目nh1和nh2来计算。

当t≤0时，无需进行移动；

当t>0时，计算平移向量vt＝[t*0.0160t*0.008]^t。

将断骨模型二按vt进行平移，即可实现粗配准，粗配准的结果如图6所示。

s4：如图7-图13对断骨模型截面部分的点集进行提取和分割：将粗配准中的特征点进行区域自生长得到扩展点集，在扩展点集中进行匹配搜索得到截面点集。本步骤将会在粗配准的基础上，对断骨模型截面部分的点集进行精确的提取和分割，分割结果如图11所示。具体步骤如下：

s41：截面特征点进行区域自生长。在粗配准过程中，可以得到最优匹配的两个聚类点集c和c′，本步骤中将这两个聚类点集反向映射回断骨模型(如图6)，可以得到两个截面的特征点集合f1和f1′，然后将这两个特征点集合进行多轮(设为5)区域自生长得到扩展后的点集f2和f2′，如图12。

s42：对f2和f2′进行匹配搜索。搜索的策略是：遍历f2中的所有数据点，计算该点处的法向量与截面处切线的夹角β，若cosβ>0.35则直接将该点提取至断面点集，否则设置一个坐标阈值r1＝20.0和向量阈值r2＝1.0，对于当前研究的数据点p，在f2′中遍历，若f2′中的点p′同时满足以下三个条件，则p和p′为匹配点：

①坐标的平方距离小于阈值，即(x1-x2)²+(y1-y2)²+(z1-z2)²<r1

②(xp+xp′)²+(yp+yp′)²+(zp+zp′)²<r2，

(xp,yp,zp)(xp′,yp′,zp′)分别是点p和p′处的法向量

③点p和p′处的法向量夹角为γ，cosγ>0。

若对点p，能找到至少5个p′满足上述条件，则将p加入至断面点集。

s43：分别使f2和f2′作为目标点集共进行两次遍历搜索，最终提取得到两个断面点集p1和p2，提取结果如图11所示。

s5：根据分割得到的两个截面点集对断骨模型进行精确配准：利用pca算法提取截面点集的主方向并根据主方向进行预配准，利用icp算法进行精配准。

预配准过程：预配准的目的是使断面大致吻合，便于后续的精配准，方法是：首先利用pca算法提取两个断面点集的主方向，以x轴增大的方向作为正方向，将两个主方向向量都调整为正方向。然后提取两个主方向向量在y-o-z平面上的投影向量v1v2，固定断骨2，将断骨1以主轴线为轴进行旋转，直至v1v2平行，此时两个断面在主轴线方向上可以做到大致吻合。预配准的结果如图13所示。

精配准过程：在预配准的基础上，利用icp算法进行精确配准，具体过程是(1)将两部分点集分别记为u与p。(2)计算最近点，即对于集合u中的每一个点，在集合p中都找出距该点最近的对应点，设集合p中由这些对应点组成的新点集为q＝{qi,i＝0,1,2,...,n}。(3)采用最小均方根法，计算点集u与q之间的配准，使得到配准变换矩阵r，t，其中r是3×3的旋转矩阵，t是3×1的平移矩阵。(4)计算坐标变换，即对于集合u，用配准变换矩阵r，t进行坐标变换，得到新的点集u1，即u1＝ru+t。(5)计算u1与q之间的均方根误差，如小于预设的极限值e，则结束，否则，以点集u1替换u，重复上述步骤。图14展示了点集精配准的效果，图15展示了两个断骨模型最终配准的结果。

s6：在断骨模型的断裂部位选取控制点，将控制点范围内的三角面片进行加厚处理得到钢板模型。用户在拼接好的断骨断面附近进行点选如图16，确定钢板模型的大概形状以及大小，记录下点选的三角平面值，选中范围内的所有表面三角面片，计算每个三角面片的法向量值，并记录下来。将每个平面根据其法向量方向进行一定程度的加厚，并填充其缝隙的部位。所得到的加厚部分即为模拟的钢板模型三维数据，可将其作为结果导出输出。图17展示了钢板拟合的效果。

本发明公开的一种弯骨骨折自动配准及内固定钢板预弯建模方法，该方法运用了曲线拟合、聚类、特征点映射等多种算法，可以高效准确地实现弯骨的自动拼接，实现了弯骨模型的全自动配准和钢板预弯，无需进行手工操作。

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，根据本发明的技术方案及其发明构思加以等同替换或改变，都应涵盖在本发明的保护范围之内。