CN111009032B - 基于改进外极线约束匹配的血管三维重建方法 - Google Patents

Abstract

本发明公开的是一种基于改进外极线约束匹配的血管三维重建方法，包括步骤：1)获取同一心动周期下两个不同视角的DSA图像的造影参数，并生成几何变换矩阵；2)分割、提取两个DSA图像中感兴趣血管段的中心线并顺序保存坐标，并对坐标进行处理来保证两中心线坐标个数相同；3)进行造影系统的自标定，优化几何变换矩阵；4)构建外极线匹配度矩阵，通过动态规划得最优点对匹配；5)由最优点对计算三维坐标，并进行三维可视化；本发明的方法对造影系统设备没有太高要求，最普遍的单面造影系统即可，且对场景和测试样本中人体拍摄角度没有太多要求，不需要附加额外的辅助设备或造作，适用于诸多场景中的血管三维重建。

Description

基于改进外极线约束匹配的血管三维重建方法

技术领域

本发明涉及一种血管三维重建方法，尤其涉及一种基于改进外极线约束匹配的血管三维重建方法，属于计算机视觉领域。

背景技术

医学图像三维重建及可视化技术能够帮助从两幅不同角度的单面造影图像重建三维冠状动脉，目前医学上主要使用冠脉造影图像技术，但是冠脉造影图像不能直接提供人体重要血管的三维空间信息，难以作到客观、准确、定量和快速的实现人体重要血管的三维重建，如采用X射线单臂造影系统对人体做X射线造影，通过旋转造影臂得到一个对应于不同造影角度的造影图像序列，单臂造影可以很方便对人体进行不同角度的造影，但是对重建来说个缺点就是我们无法得到同一时刻的不同视角的造影图像，这给重建带来很大的困难。要重建血管的真实三维空间结构，需要得到血管至少两个不同角度的投影信息，传统方法首先提取出血管的骨架，然后通过不同视角空间约束关系，对不同视角投影图像的血管像素点通过立体视觉中的外极线约束进行正确匹配并重建，但当血管比较扭曲或外极线不精确时，外极线会与血管有多个交点。

目前常见的寻找共轭点方法主要有采用多种匹配度评价的混合匹配策略，如公开号CN104361626A的中国专利申请文献通过单应性矩阵及外极线约束的混合匹配重建皮下静脉，具体为：首先通过外极线约束构建候补匹配点集，再通过寻找SURF匹配特征点后计算两视图问的单应矩阵，从匹配点集中寻找匹配差异度小于匹配阈值的点。然而对于心血管这种空间立体性较强的血管并不适用单应矩阵，又如公开号CN101105393A的中国专利申请文件公开的一种基于投射光栅的物体表面三维轮廓重建后测量，包括以下步骤：通过向物体表面投射不同频率的光栅，获得多个相位图；采用频率合成方法，展开相位图；利用展开的相位与外极线匹配对应点，计算物体表面的三维轮廓。这是人为添加了额外的信息，对于造影图像中剩余可用的信息只剩血管半径，但血管半径在血管中变化变化不大，且由于血管横截面是不规则图形，因此同一血管点在不同视角下的半径也可能不同，所以造影图像中并没有真正可行的额外信息来辅助外极线约束进行混合匹配。

现阶段基于外极线约束匹配的血管三维重建方法，要求选择相同心动周期的图像，且对于外极线的计算精度要求较高，对于出现外极线与血管有多个交点时，可能无法判断出其中真正的匹配点。

发明内容

为了解决现有技术中存在的问题，本发明提供一种基于改进外极线约束匹配的血管三维重建方法，可以在不用人为添加额外信息的情况下，解决外极线匹配时外极线与血管有多个交点的问题，较好地提升外极线匹配的精度，从而提升血管三维重建的精度。

为了实现上述目的，本发明是通过以下技术方案实现的：

基于改进外极线约束匹配的血管三维重建方法，该血管三维重建方法包括以下步骤：

1)读取同一个人在不同角度下的造影序列，记录冠脉造影系统的成像的相关参数，获取相应心动周期下同一心动时刻的DSA图像及造影参数，并生成几何变换矩阵GT；针对目前临床上普遍采用单面造影系统的情况，本发明主要适用于单面冠脉造影图像；

2)获取两个DSA图像中感兴趣血管段的中心线坐标，并对坐标进行处理，保证两中心线坐标个数相同；实际两条中心线都可以保存为连通的、有序的点序列，但由于中心线长度不同，会在之后的外极线匹配对应点时产生误差，因此通过插值的方法对中心线进行插值，使得两条中心线的长度相等，从而原则上能做到中心线上的点一一对应；

3)进行造影系统的自标定，构建目标函数，优化几何变换矩阵，确保生成精准的外极线；所述造影系统自标定是利用造影图像自身的特点，无需附加模型或道具，也不改变传统造影程序，简单易行，可以达到较高的标定精度；

4)计算外极线匹配度，通过动态规划获得最优点对匹配；

5)由最优点对计算三维坐标并显示，三维坐标生成方法采用计算速度更快的解析方法，最小二乘法。

作为一种优选，所述血管三维重建方法适用于冠状动脉血管以及其他血管的三维重建，所述其他血管包括肾动脉、股动脉、下腔静脉。

作为一种优选，在步骤1)中，通过造影时同步记录的心电信号来保证同一心动周期，根据心电信号选择对应心动周期内同一时刻的造影图像，由于在心脏舒张最大时心脏运动速率相对比较小，局部的复杂运动较少，更容易观察冠状动脉的结构和形态，所以一般选择心脏舒张最大时刻的造影图像用于冠状动脉的三维重建；

作为一种优选，步骤1)中所述造影参数包括内参、外参，所述内参包括：像素间距q，成像大小W×H，射线源到投影面中心距离SID；所述外参包括：射线源到同心点距离SOD，造影左右角度LAO/RAO，造影前后角度CRAN/CAUD。

作为一种优选，步骤1)中几何变换矩阵GT包括R＝R_x(β₂)·R_y(α₂)·R_y(-α₁)·R_x(-β₁)， 其中，所述/>表示绕Y轴顺时针旋转θ角度的旋转矩阵，所述/>表示绕X轴顺时针旋转θ角度的旋转矩阵，所述/>表示从射线源0点沿Z轴平移L的平移矩阵。

作为一种优选，步骤2)中所述中心线坐标以首尾顺序存储。

作为一种优选，比较不同造影角度的中心线上投影点个数，以多的投影点个数作为标准，将投影点少的中心线用插值成相同个数的投影点。

作为一种优选，所述目标函数为图像二维重建误差，具体为图像投影点与三维重建后反投影点之间的欧式距离，所用的优化方法为最优化理论中的相关方法，所述相关方法包括最速下降法、Levenberg-Marquard。

作为一种优选，步骤4)中通过动态规划获得最优点对匹配，由于在应用外极线匹配时，因为血管的扭曲或者视角的变换，往往会造成一条外极线与，血管中心线有多个交点，使得某一点有多个对应的匹配点，本发明通过构建全局的外极线误差矩阵，考虑每一步前进的最优性，从而得到全局最优匹配路径，具体过程包括：

a)构建候补匹配点集，通过外极线约束，对候补匹配点集中的所有点对计算外极线匹配度Degree；

b)将匹配度Degree依次存入n*n矩阵中，构成匹配度矩阵，其中n为中心线投影点个数；

c)用动态规划的方法求出从匹配度矩阵左上角到矩阵右下角的最短路径，该动态规划的策略是从边界开始，逐段递推寻优，直到走到另一个边界，从候补匹配点集中筛选出最优匹配点对。

有益效果：本发明的方法对造影系统设备要求低，最普遍的单面造影系统即可；对场景和测试样本中人体拍摄角度没有太多要求，不需要附加额外的辅助设备或造作，适用于诸多场景中的血管三维重建；可以在不用人为添加额外信息的情况下，解决外极线匹配时外极线与血管有多个交点的问题，较好地提升外极线匹配的精度，从而提升血管三维重建的精度。

附图说明

图1为血管重建方法流程图。

图2为造影图像一。

图3为造影图像二。

图4为造影图像一中心线提取示意图。

图5为造影图像二中心线提取示意图。

图6为造影图像一外极线匹配示意图。

图7为造影图像二外极线匹配示意图。

图8为三维重建模型第三种角度展示图。

图9为三维重建模型第四种角度展示图。

具体实施方式

以下结合说明书附图，对本发明作进一步说明，但本发明并不局限于以下实施例。

具体的：下面结合实施例和附图1-9来详细说明本发明，并将本发明的实验中相关数据公开阐述，进而将本发明的方法详细透彻的传达给本来领域的技术人员。

本发明针对利用外极线约束获得血管对应点匹配时有多个交点的问题，采用矩阵描述外极线误差，通过动态规划的方法最优化匹配路径，从而可以精确的进行三维重建。

如图1所示，本实施例中的具体实现步骤如下：

1)读取同一个人在不同角度下的造影序列，并记录冠脉造影系统的成像的相关参数，读取对应DICOM文件中的ECG数据，分析其心电信号，选择对应心动周期内心脏舒张最大时刻的造影图像，即在心电信号为波峰时记录的造影图像，所述造影系统如图2所示，采用单面冠脉造影系统；

2)在选定的造影图像一(图2所示)、造影图像二(图3所示)中根据血管树的拓扑结构，选择对应的分叉点分别作为感兴趣血管的起始点和终止点，运用现有技术成熟的基于区域生长的图像分割算法和Zhang Suen细化算法，对血管进行分割、细化后提取血管中心线，然后顺序存储中心线点坐标，如图4、5所示，最后对存储的血管中心线进行插值处理。具体为统计血管中心线一(图4所示)、中心线二(图5所示)上点的数量，选取中心线进行插值，保证两中心线上点的数量一致；

3)对造影系统进行自标定，具体包括如下步骤：

a)首先通过记录的相关参数构建内参矩阵及初始几何变换矩阵GT；

造影系统中的内部参数：像素间距q，成像大小W×H，射线源到投影面中心距离SID，用于将X射线源局部坐标系中三维点(x_i，y_i，z_i)投影到图像坐标系上的投影点(e_i，f_i)，转换关系如下：

造影系统中的外部参数：射线源到同心点距离SOD，造影左右角度LAO/RAO，造影前后角度CRAN/CAUD。

几何变换矩阵GT中R＝R_x(β₂)·R_y(α₂)·R_y(-α₁)·R_x(-β₁)，其中α₁为造影图像一中LAO/RAO角度，RAO为正，β₁为造影图像一中CRAN/CAUD角度，CRAN为正，α₂、β₂为造影图像二对应的角度；

中的L₁为造影图像一中的SOD，/>中的L₂为造影图像二中的SOD；

b)选取中心线一、中心线二中各自的起始点和终止点作为数据源，选取合适的误差作为目标函数，以最速下降法或Levenberg-Marquard等作为优化方法，得到优化后的内参矩阵及几何变换矩阵GT，即完成了造影系统的自标定；

4)计算外极线匹配度，构建n*n的匹配度矩阵；

外极线的计算过程如下，对于造影图像一上的点p₁(u₁，v₁)在造影图像二上对应的外极线L₂方程可表示为：

ξ₂·(a₃·b₂-a₂·b₃)+η₂·(a₁·b₃-a₃·b₁)+(a₂·b₁-a₁·b₂)＝0；即造影图像二上的对应点q₁(u₂，v₂)满足L₂方程；

5)由于在应用外极线匹配时，面对血管中心线与外极线有多个交点的情况，如图6、7所示，对构建的匹配度矩阵进行全局的动态规划，获得从起始点到终止点的最优匹配路径；

6)最后通过前面优化后的最优匹配点对计算三维坐标，本发明采用最小二乘法计算：

有方程

写成A·C＝B，该方程由四个线性方程组组成，求解3个未知量x₁，y₁，z₁，因此是一个超限定方程组，可通过最小二乘法求解。即若已知造影图形一、二的对应点对，可计算出该三维点在X射线源一坐标系X₁Y₁Z₁S₁中的坐标(x₁，y₁，z₁，采用VTK工具包对三维点进行三维可视化，如图8、9所示，图8为三维重建模型第三种角度展示图，图9为三维重建模型第四种角度展示图。

最后，需要注意的是，本发明不限于以上实施例，还可以有很多变形。本领域的普通技术人员能从本发明公开的内容中直接导出或联想到的所有变形，均应认为是本发明的保护范围。

Claims (8)

1.基于改进外极线约束匹配的血管三维重建方法，其特征在于该血管三维重建方法包括以下步骤：

1)读取同一个人在不同角度下的造影序列，获取相应心动周期下同一心动时刻的DSA图像及造影参数，并生成几何变换矩阵GT；记录冠脉造影系统的成像的相关参数，读取对应DICOM文件中的ECG数据，分析其心电信号，选择对应心动周期内心脏舒张最大时刻的造影图像，即在心电信号为波峰时记录的造影图像；

2)获取两个DSA图像中感兴趣血管段的中心线坐标，并对坐标进行处理，保证两中心线坐标个数相同；在选定的造影图像一、造影图像二中根据血管树的拓扑结构，选择对应的分叉点分别作为感兴趣血管的起始点和终止点，运用基于区域生长的图像分割算法和ZhangSuen细化算法，对血管进行分割、细化后提取血管中心线，然后顺序存储中心线点坐标，最后对存储的血管中心线进行插值处理；具体为统计血管中心线一、中心线二上点的数量，选取中心线进行插值，保证两中心线上点的数量一致；

3)进行造影系统的自标定：构建目标函数，优化几何变换矩阵，生成精准的外极线；a)首先通过记录的相关参数构建内参矩阵及初始几何变换矩阵GT；

造影系统中的内部参数：像素间距q，成像大小W×H，射线源到投影面中心距离SID，用于将X射线源局部坐标系中三维点(x_i，y_i，z_i)投影到图像坐标系上的投影点(e_i，f_i)，转换关系如下：

造影系统中的外部参数：射线源到同心点距离SOD，造影左右角度LAO/RAO，造影前后角度CRAN/CAUD；

几何变换矩阵GT中R＝R_x(β₂)·R_y(α₂)·R_y(-α₁)·R_x(-β₁)，其中α₁为造影图像一中LAO/RAO角度，RAO为正，β₁为造影图像一中CRAN/CAUD角度，CRAN为正，α₂、β₂为造影图像二对应的角度；

中的L₁为造影图像一中的SOD，/>中的L₂为造影图像二中的SOD；

b)选取中心线一、中心线二中各自的起始点和终止点作为数据源，选取合适的误差作为目标函数，以最速下降法或Levenberg-Marquard作为优化方法，得到优化后的内参矩阵及几何变换矩阵GT，即完成了造影系统的自标定；

4)计算外极线匹配度，通过动态规划获得最优点对匹配；

计算外极线匹配度，构建n*n的匹配度矩阵；

外极线的计算过程如下，对于造影图像一上的点p₁(u₁，v₁)在造影图像二上对应的外极线L₂方程可表示为：

ξ₂·(a₃·b₂-a₂·b₃)+η₂·(a₁·b₃-a₃·b₁)+(a₂·b₁-a₁·b₂)＝0；

造影图像二上的对应点q₁(u₂，v₂)满足L₂方程；

5)由最优点对计算三维坐标并显示；

由于在应用外极线匹配时，面对血管中心线与外极线有多个交点的情况，对构建的匹配度矩阵进行全局的动态规划，获得从起始点到终止点的最优匹配路径；最后通过前面优化后的最优匹配点对计算三维坐标，采用最小二乘法计算：

有方程

写成A·C＝B，该方程由四个线性方程组组成，计算出该三维点在X射线源一坐标系X₁Y₁Z₁S₁中的坐标(x₁，y₁，z₁)，采用VTK工具包对三维点进行三维可视化。

2.根据权利要求1所述的基于改进外极线约束匹配的血管三维重建方法，其特征在于：所述血管三维重建方法适用于冠状动脉血管以及其他血管的三维重建，所述其他血管包括肾动脉、股动脉、下腔静脉。

3.根据权利要求1或2所述的基于改进外极线约束匹配的血管三维重建方法，其特征在于：在步骤1)中，通过造影时同步记录的心电信号来保证同一心动周期。

4.根据权利要求1所述的基于改进外极线约束匹配的血管三维重建方法，其特征在于：步骤1)中几何变换矩阵GT包括R＝R_x(β₂)·R_y(α₂)·R_y(-α₁)·R_x(-β₁)， 其中，所述表示绕Y轴顺时针旋转θ角度的旋转矩阵，所述表示绕X轴顺时针旋转θ角度的旋转矩阵，所述/>表示从射线源O点沿Z轴平移L的平移矩阵。

5.根据权利要求1或2或4所述的基于改进外极线约束匹配的血管三维重建方法，其特征在于：步骤2)中所述中心线坐标以首尾顺序存储。

6.根据权利要求3所述的基于改进外极线约束匹配的血管三维重建方法，其特征在于：比较不同造影角度的中心线上投影点个数，以多的投影点个数作为标准，将投影点少的中心线用插值成相同个数的投影点。

7.根据权利要求1所述的基于改进外极线约束匹配的血管三维重建方法，其特征在于：所述目标函数为图像二维重建误差，具体为图像投影点与三维重建后反投影点之间的欧式距离，所用的优化方法为最优化理论中的相关方法，所述相关方法包括最速下降法、Levenberg-Marquard。

8.根据权利要求1所述的基于改进外极线约束匹配的血管三维重建方法，其特征在于：步骤4)中通过动态规划获得最优点对匹配，具体过程包括：

a)构建候补匹配点集，通过外极线约束，对候补匹配点集中的所有点对计算外极线匹配度Degree；

b)将匹配度Degree依次存入n*n矩阵中，构成匹配度矩阵，其中n为中心线投影点个数；

c)用动态规划的方法求出从匹配度矩阵左上角到矩阵右下角的最短路径，该动态规划的策略是从边界开始，逐段递推寻优，直到走到另一个边界，从候补匹配点集中筛选出最优匹配点对。