CN107316327A - 基于最大公共子图与包围盒的断骨断面及断骨模型配准方法 - Google Patents

Abstract

本发明公开了一种基于最大公共子图与包围盒的断骨断面及断骨模型配准方法，包括如下步骤：断骨轴线提取步骤：断骨高斯映射步骤：断面点集提取步骤：通过剔除其法向量与所述断骨轴线呈一定角度内的三角面片，将所述的断骨模型分割成两个断骨模型；每个模型包含该断骨的断面模型；分别对两个断骨进行三角面片的剔除，最终得到四个子模型，其中2个为断面模型；提取4个所述子模型包含节点数的断面点集，通过使用最大公共子图算法两两比较所述的断面电机，选择2个包含节点数最大的一对子模型作为所述的两个断面模型；断骨预配准步骤：建立所述的两个断面模型的包围盒，将两包围盒在空间上对齐完成预配准；断骨精确配准步骤：通过ICP算法对两断骨模型进行精确配准，完成两断骨的接合。

Description

基于最大公共子图与包围盒的断骨断面及断骨模型配准方法

技术领域

本发明涉及一种基于最大公共子图与包围盒的断骨断面及断骨模型配准方法。涉及专利分类号G06计算；推算；计数G06T一般的图像数据处理或产生 G06T7/00图像分析，例如从位像到非位像。

背景技术

目前骨折手术普遍采用人工复位和伤肢内固定相结合的方法，这种方法存在的问题是创伤大、出血多、容易引发神经血管损伤等并发症。因此，可以利用计算机对断骨模型进行虚拟拼接，从而在术前得到钢板的各种几何参数。然而目前的断骨虚拟拼接方法存在很多问题，比如断面分割不准确、需要手动操作、预配准不精确等等，这些问题都很大程度地影响了断骨模型配准的效果。

在现有技术中，存在通过根据断面法向量与主轴线夹角和法向量的突变两个方面来进行断面分割的手段，也存在直接根据断骨的轴线进行预配准的方案。但是，根据断面法向量与主轴线夹角和法向量的突变两个方面来进行断面分割，需要人为的在模型上进行点选，不能实现断面的自动提取。

另外，单纯将断骨的主轴线进行对齐并不能满足预配准的要求，可能会出现断骨的截面方向相背离的情况，而且截面形状也无法做到大致吻合，很容易造成精配准的不准确。

发明内容

本发明要解决的技术问题是：一种基于最大公共子图与包围盒的断骨断面及断骨模型配准方法，包括如下步骤：

断骨轴线提取步骤：

根据经由CT图像序列重建生成的断骨模型，生成该模型初始样本矩阵Z；对该样本矩阵进行中心化，得到矩阵X；计算该矩阵X的协方差矩阵C和协方差矩阵C的特征值，选取最大特征值对应的特征向量作为断骨轴线；

断骨高斯映射步骤：

计算所述断骨模型上每一个三角面片的法向量，对全部的法向量进行单位化，得到每一个三角面片的单位法向量n；将得到的所有单位法向量n的起点平移到坐标原点，向量的终点则会落到以坐标原点O为球心，半径为1的单位球面S上，断骨模型的高斯映射完成；

断面点集提取步骤：

通过剔除其法向量与所述断骨轴线呈一定角度内的三角面片，将所述的断骨模型分割成两个断骨模型；每个模型包含该断骨的断面模型；

分别对两个断骨进行三角面片的剔除，最终得到四个子模型，其中2个为断面模型；

提取4个所述子模型包含节点数的断面点集，通过使用最大公共子图算法两两比较所述的断面电机，选择2个包含节点数最大的一对子模型作为所述的两个断面模型；

断骨预配准步骤：

建立所述的两个断面模型的包围盒，将两包围盒在空间上对齐完成预配准；

断骨精确配准步骤：

通过ICP算法对两断骨模型进行精确配准，完成两断骨的接合。

作为优选的实施方式，所述的断骨轴线提取步骤具体包括如下步骤：

a.对断骨模型中的n个顶点，将所述的n个点坐标赋值给一个3×n的矩阵 Z，生成该模型的初始样本矩阵 为三维列向量；

b.对所述的初始样本矩阵进行中心化，得到中心化后的矩阵X；

其中：计算得到的三维向量则可以视为断骨模型在空间上的中心；

c.计算所述矩阵的协方差阵列C，转换过程如下公式所示：

计算所述的协方差阵列C为3×3矩阵；

d.分别求出所述协方差阵列C的特征值λ₁，λ₂，λ₃和特征值对应的特征向量

比较所述的特征值λ₁，λ₂，λ₃，选择最大的特征值对应的特征向量，作为该断骨模型的轴线向量。

作为优选的实施方式，所述的断骨模型的高斯映射步骤具体包括如下步骤：

a.计算所述三维模型上每个三角面片的法向量；

b.对得到的所述每个三角面片的法向量进行单位化，得到每个所述三角面片的单位法向量n；

c.将得到的每个单位方向量的起点平移至坐标原点，则每个所述的单位法向量的重点最终落在以坐标原点O为球心，半径为1的单位球面S上，完成高斯映射。

作为优选的实施方式，所述的断面点集提取步骤中剔除的过程如下：

计算每个所述三角面片法向量与所述断骨模型的轴线方向向量的夹角；剔除法向量所成夹角在0.4π-0.6π之间的三角面片，将剩余的法向量与轴线方向向量夹角范围在[0,0.4π]∪[0.6π,π]内的三角面片保留，将一个断骨模型分割成2 个子模型，其中一个所述的子模型为断面模型；

对所述的2个断骨模型中将夹角在0.4π-0.6π之间的三角面片剔除，完成断骨模型侧面的三角面片剔除，得到四个子模型，其中两个是断面模型；

更进一步的，得到两个断面模型的过程如下：

使用最大公共子图算法判断断面模型，分别对一个断骨模型的两个子模型与另一断骨的两个子模型使用最大公共子图算法，最大公共子图所包含节点数最大的一对子模型，即为2个断面模型。

更进一步的，所述的精确配准断骨精确配准步骤过程如下：

将两断面模型分点集分别记为P与Q；在目标点集P中选取点集合P_i ^k∈P，计算源点集合Q中相应的点满足

计算旋转矩阵R^k和平移矢量T^k，满足计算 P^k+1＝{P_i ^{k +1}|P_i ^k+1＝R^kP_i ^k+T^k,P_i ^k∈P}和其中，K代表着第k次迭代，P_i ^k和分别代表目标点集和源点集中的点，min代表设定好的最小距离的阈值，d^k+1代表着第k+1次迭代下的平均最小距离；

若d^k+1大于等于事先设定的阈值t，则重新在目标点集P中选取点集合步骤，重新开始算法，当d^k+1<t或循环次数大于事先设定好的循环次数的阈值时，跳出循环。

作为优选的实施方式，断骨精确配准步骤之后还具有虚拟钢板预弯步骤：

在拼接好的断骨断面附近进行点选，确定钢板模型的形状以及尺寸；

记录下点选的三角平面值，选中范围内的所有表面三角面片；

计算每个三角面片的法向量值；将每个平面根据其法向量方向进行一定程度的加厚，并填充其缝隙的部位；

所得到的加厚部分即为模拟的钢板模型三维数据。

附图说明

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

图1为本发明的算法流程图

图2为本发明断骨轴线提取后的效果示意图

图3为本发明断骨模型高斯映射的结果示意图

图4为本发明断骨模型分割成4个子模型的分割结果示意图

图5为本发明使用最大公共子图方法提取的断面模型示意图

图6为本发明基于包围盒对齐的断骨模型预配准结果示意图

图7为本发明断骨模型精确配准结果示意图

图8为本发明断骨表面选取的特征点示意图

图9为本发明虚拟钢板模型示意图

具体实施方式

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

如图1-9所示：基于最大公共子图与包围盒的断骨断面及断骨模型配准方法，主要包括如下步骤：

断骨轴线提取

本方法利用PCA算法进行断骨的轴线提取。提取效果如图2所示。

输入初始样本矩阵。根据输入的模型，设模型一共包含n个点，将这n个点的坐标值赋值给一个3×n的矩阵，此矩阵则为初始样本矩阵

初始样本矩阵中心化。将初始样本矩阵Z中心化后，得到矩阵中心化的公式为其中，计算得到的三维向量则可以视为断骨模型在空间上的中心。

计算矩阵X的协方差阵列C。C的公式为计算得到的C是一个 3×3的矩阵。

求协方差矩阵C的特征值和特征向量。对于得到的协方差阵列C，求出其特征值λ₁，λ₂，λ₃与对应的特征向量将三个特征值做比较，得到最大的特征值，其对应的特征向量，就是该断骨模型的轴线方向。

断骨高斯映射

本方法将断骨模型进行了高斯映射操作，映射结果如图3所示。计算断骨模型上的每一个三角面片的法向量。

将计算得到的三角面片的法向量单位化，单位化的过程可采用如下所述的方法：

设法向量为e＝(e1，e2，e3)，那个其单位法向量

得到每个三角面片的单位法向量n。将得到的所有单位法向量的起点平移到坐标原点，向量的终点则会落到以坐标原点O为球心，半径为1的单位球面S 上，断骨模型的高斯映射完成。

通过在高斯映射中将每个面片的法向量计算出来，通过对高斯映射的显示，可以直观的观察到法向量与模型轴线的关系，为下一步根据法向量与轴线夹角的关系剔除三角面片做铺垫。

断骨断面数据点集的提取

断骨模型分割。根据每个三角面片的法向量，计算其与模型轴线方向向量的夹角，将与轴线方向向量夹角在(0.4π,0.6π)范围内的三角面片剔除，则留下的三角面片法向量与轴线方向向量的夹角范围是[0,0.4π]∪[0.6π,π]，这样，就把一个断骨模型分割成了两个子模型，其中包含该断骨的断面模型。分别对两个断骨进行三角面片的剔除，则可以获得四个子模型，其中有两个是断面模型。模型分割结果如图4所示。

基于最大公共子图的断面数据点集提取。本方法提出了一种基于最大公共子图算法的方法来提取断面点集。提取结果如图5所示。本方法的大致实现步骤是：

构建图结构。将模型的一个三角面片看作是图结构中的一个顶点，如果两个三角面片在模型中是边相邻的关系，那么代表两个三角面片的顶点之间则存在一条边将这两个顶点连接起来。用该顶点所表示的三角面片的法向量与其边相邻的三个三角面片的法向量的夹角之和来表示每个顶点的值。对于每条边的值，只要两个顶点之间有边的关系，则将该边的值赋值为1。

使用最大公共子图算法判断断面模型。分别对每个子模型与另一个模型的两个子模型使用最大公共子图算法，最大公共子图所包含的节点数最大的一对子模型，即为两个断面模型。

最大公共子图算法流程如下：

(1)对于已经建立好的两个图G₁、G₂，构建关联图G；

(2)根据图G的结构，分别对remaining\connected和remaining∩connected两个节点集进行染色，得到colourClasses；

(3)对colourClasses中每个颜色类中的每一个节点v进行判断，如果当前的结果集solution的节点个数与colourClasses的类的个数小于最优解incumbent的节点个数或者当前节点不在connected点集里并且solution集合不为空集，则删除当前颜色类，从下一个颜色类进行节点遍历；

(4)否则，将当前节点v添加进solution集合里，如果新的solution的节点个数大于目前最优解incumbent的节点个数，则更新incumbent集合，并更新connected 和remaining集合；

(5)如果当前remaining集合不为空，则向下递归寻找最大公共子图，直到remaining集合为空集。

对于以上流程，其中，solution集合中包含当前已经寻找到的公共子图的节点，remaining集合中包含与solution中的节点相邻的和还没有被接受或者拒绝的节点，connected集合中包含在G₁中与已经被接受的顶点相邻的顶点所匹配的节点对。

基于包围盒对齐的断骨预配准

断骨模型的预配准，是通过将两个断面模型的包围盒在空间上对齐的方式来实现的。预配准结果如图6所示。

将包围盒的短轴方向固定在Y轴上，长轴方向固定在X轴上，中轴方向固定在Z轴方向。如果发现某一个包围盒的长轴方向与中轴方向分别在Z轴和X 轴上，则需要对该包围盒以及对应的断骨模型进行空间变换。先是将包围盒以及断骨模型绕Y轴旋转90度，再沿着X轴方向平移包围盒长轴的长度，这样就可以实现使包围盒三个轴的方向分别于三个坐标轴对应。

经过第一步变换之后，断骨模型会有两种状态，一种是整个断骨模型完全在XOZ平面的一侧，另一种是断骨模型被XOZ平面分为两部分。对于第二种情况，我们还需要对模型位置进行处理，将其平移到XOZ平面的一侧。对于平移的方向与距离，如果断骨模型的中心点的Y坐标大于0，则将模型沿着Y轴正方向平移包围盒短轴长度的距离，如果断骨模型的中心点的Y坐标小于0，则将模型沿Y轴负方向平移包围盒短轴长度的距离。

两个断骨模型被变换到了XOZ平面的同侧或异侧。对于在XOZ平面同侧的两个断骨模型，需要将其变换到分别在XOZ平面的两侧。变换的过程为，将一个断骨模型绕着X轴或者Z轴旋转180度，再沿着Z轴方向或者X轴方向平移中轴或者长轴长度的距离。

将两个断骨模型的中心点坐标与包围盒的中心点坐标进行比较，如果包围盒中心点的X坐标在两个断骨模型中心点坐标的中间，则两个断骨已经预配准完成；如果包围盒中心点的X坐标比两个断骨模型中心点的X坐标都大或者都小，则说明发生了180度翻转的情况。此时，只需要将一个模型绕着Y轴旋转 180度，并在X轴方向平移包围盒长轴长度，在Z轴方向平移包围盒中轴长度即可。

断骨精配准

预配准结束后，利用ICP算法对两组断骨模型点集进行精配准操作，进一步提高其准确性。精配准结果如图7所示。

将两部分点集分别记为P与Q。

在目标点集P中选取点集合P_i ^k∈P，计算源点集合Q中相应的点满足

计算旋转矩阵R^k和平移矢量T^k，满足

计算P^k+1＝{P_i ^k+1|P_i ^k+1＝R^kP_i ^k+T^k,P_i ^k∈P}和

若是d^k+1大于等于事先设定的阈值t，则跳转到步骤2，当d^k+1<t或循环次数大于事先设定好的循环次数的阈值时，跳出循环。

虚拟钢板预弯

在本步骤，主要利用对断骨断裂部位曲面数据的相应操作来模拟出钢板的大概形态。

在拼接好的断骨断面附近进行点选(如图8)，确定钢板模型的大概形状以及大小。

记录下点选的三角平面值，选中范围内的所有表面三角面片。计算每个三角面片的法向量值，并记录下来。将每个平面根据其法向量方向进行一定程度的加厚，并填充其缝隙的部位。所得到的加厚部分即为模拟的钢板模型三维数据(如图9)，可将其作为结果导出输出。

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

Claims (7)

1.一种基于最大公共子图与包围盒的断骨断面及断骨模型配准方法，其特征在于包括如下步骤：

断骨轴线提取步骤：

根据经由CT图像序列重建生成的断骨模型，生成该模型初始样本矩阵Z；对该样本矩阵进行中心化，得到矩阵X；计算该矩阵X的协方差矩阵C和协方差矩阵C的特征值，选取最大特征值对应的特征向量作为断骨轴线；

断骨高斯映射步骤：

计算所述断骨模型上每一个三角面片的法向量，对全部的法向量进行单位化，得到每一个三角面片的单位法向量n；将得到的所有单位法向量n的起点平移到坐标原点，向量的终点则会落到以坐标原点O为球心，半径为1的单位球面S上，断骨模型的高斯映射完成；

断面点集提取步骤：

通过剔除其法向量与所述断骨轴线呈一定角度内的三角面片，将所述的断骨模型分割成两个断骨模型；每个模型包含该断骨的断面模型；

分别对两个断骨进行三角面片的剔除，最终得到四个子模型，其中2个为断面模型；

提取4个所述子模型包含节点数的断面点集，通过使用最大公共子图算法两两比较所述的断面点集，选择2个包含节点数最大的一对子模型作为所述的两个断面模型；

断骨预配准步骤：

建立所述的两个断面模型的包围盒，将两包围盒在空间上对齐完成预配准；

断骨精确配准步骤：

通过ICP算法对两断骨模型进行精确配准，完成两断骨的接合。

2.根据权利要求1所述的基于最大公共子图与包围盒的断骨断面及断骨模型配准方法，其特征还在于：所述的断骨轴线提取步骤具体包括如下步骤：

a.对断骨模型中的n个顶点，将所述的n个点坐标赋值给一个3×n的矩阵Z，生成该模型的初始样本矩阵 为三维列向量；

b.对所述的初始样本矩阵进行中心化，得到中心化后的矩阵X；

其中：计算得到的三维向量则可以视为断骨模型在空间上的中心；

c.计算所述矩阵的协方差阵列C，转换过程如下公式所示：

计算所述的协方差阵列C为3×3矩阵；

d.分别求出所述协方差阵列C的特征值λ₁，λ₂，λ₃和特征值对应的特征向量

比较所述的特征值λ₁，λ₂，λ₃，选择最大的特征值对应的特征向量，作为该断骨模型的轴线向量。

3.根据权利要求1所述的基于最大公共子图与包围盒的断骨断面及断骨模型配准方法，其特征还在于所述的断骨模型的高斯映射步骤具体包括如下步骤：

a.计算所述三维模型上每个三角面片的法向量；

b.对得到的所述每个三角面片的法向量进行单位化，得到每个所述三角面片的单位法向量n；

c.将得到的每个单位方向量的起点平移至坐标原点，则每个所述的单位法向量的重点最终落在以坐标原点O为球心，半径为1的单位球面S上，完成高斯映射。

4.根据权利要求1所述的基于最大公共子图与包围盒的断骨断面及断骨模型配准方法，其特征还在于所述的断面点集提取步骤中剔除的过程如下：

计算每个所述三角面片法向量与所述断骨模型的轴线方向向量的夹角；

剔除法向量所成夹角在0.4π-0.6π之间的三角面片，将剩余的法向量与轴线方向向量夹角范围在[0,0.4π]∪[0.6π,π]内的三角面片保留，将一个断骨模型分割成2个子模型，其中一个所述的子模型为断面模型；

对所述的2个断骨模型中将夹角在0.4π-0.6π之间的三角面片剔除，完成断骨模型侧面的三角面片剔除，得到四个子模型，其中两个是断面模型。

5.根据权利要求4所述的基于最大公共子图与包围盒的断骨断面及断骨模型配准方法，其特征还在于得到两个断面模型的过程如下：

使用最大公共子图算法判断断面模型，分别对一个断骨模型的两个子模型与另一断骨的两个子模型使用最大公共子图算法，最大公共子图所包含节点数最大的一对子模型，即为2个断面模型。

6.根据权利要求5所述的基于最大公共子图与包围盒的断骨断面及断骨模型配准方法，其特征还在于所述的精确配准断骨精确配准步骤过程如下：

将两断面模型分点集分别记为P与Q；在目标点集P中选取点集合P_i ^k∈P，计算源点集合Q中相应的点满足

计算旋转矩阵R^k和平移矢量T^k，满足计算P^k+1＝{P_i ^k+1|P_i ^k+1＝R^kP_i ^k+T^k,P_i ^k∈P}和其中，K代表着第k次迭代，P_i ^k和分别代表目标点集和源点集中的点，min代表设定好的最小距离的阈值，d^k+1代表着第k+1次迭代下的平均最小距离；

若d^k+1大于等于事先设定的阈值t，则重新在目标点集P中选取点集合步骤，重新开始算法，当d^k+1<t或循环次数大于事先设定好的循环次数的阈值时，跳出循环。

7.根据权利要求1所述的基于最大公共子图与包围盒的断骨断面及断骨模型配准方法，其特征还在于断骨精确配准步骤之后还具有虚拟钢板预弯步骤：

在拼接好的断骨断面附近进行点选，确定钢板模型的形状以及尺寸；

记录下点选的三角平面值，选中范围内的所有表面三角面片；

计算每个三角面片的法向量值；将每个平面根据其法向量方向进行一定程度的加厚，并填充其缝隙的部位；

所得到的加厚部分即为模拟的钢板模型三维数据。