CN111260704B - 基于启发式树搜索的血管结构3d/2d刚性配准方法及装置 - Google Patents

Abstract

一种基于启发式树搜索的血管结构3D/2D刚性配准方法及装置，其3D图像和2D图像的血管图匹配精确度高。方法包括：(1)利用3D和2D血管拓扑一致性来实现血管图匹配，将图的匹配结果表示为成对的3D和2D边的集合；(2)将血管匹配过程看作是在已有的血管对上增加一对新的匹配边的连续过程，将匹配过程分解为连续状态，并用于搜索树的构造；(3)在搜索树的每个节点上，使用一个封闭解来计算基于血管点密集匹配的配准结果，设计一个评价配准质量和配准质量的节点评分；(4)在A‑star搜索算法的基础上，利用一种改进的启发式树搜索策略，寻找节点得分最高的最优结果。

Description

基于启发式树搜索的血管结构3D/2D刚性配准方法及装置

技术领域

本发明涉及医学图像处理的技术领域，尤其涉及一种基于启发式树搜索的血管结构3D/2D刚性配准方法，以及基于启发式树搜索的血管结构3D/2D刚性配准装置。

背景技术

目前，微创介入治疗是血管疾病的主要治疗手段，在介入治疗过程中，手术器械的操作以X射线血管造影(XRA)图像为指导。造影剂通过导管注入感兴趣的动脉并进行成像，这种成像模式在显示血管管腔方面有令人满意的表现。在XRA中，血管内导航的手术器械也可以清楚地显示出来。然而，由于XRA缺乏空间信息，在单视图二维投影的指导下进行准确的介入操作对介入医师来说是困难的。因此，在介入手术中，医生经常使用旋转C型臂获得的多视图血管造影图像，然而这样会加大造影剂的注入且会给患者带来负担。为了解决这个问题，可以将术前计算机断层扫描血管造影(CTA)图像与术中XRA成像相结合使用。通过将三维血管投影覆盖在二维实时图像上来增强介入图像，医生和病人均可以从多模态数据融合和可视化中受益。为了实现这一目标，3D/2D配准技术是其中获得良好的对齐和对应关系的关键。

3D/2D配准方法通常使用一张术前的3D图像和作为配准源的多张术中2D X射线造影图来实现配准。3D和2D图像数据以及C型臂的成像几何参数需要作为配准模型的输入。对于X射线图像序列，通常认为3D/2D配准是3D/2D+t的。旋转成像是临床常用的一种多平面X射线成像技术，在这种情况下，多平面3D/2D配准可视为3D/2.5D。由于3D/2D配准是3D/2D+t和3D/2.5D的基础，因此本文只讨论3D图像到单帧单平面2D图像的配准。根据配准技术的性质，将3D/2D配准方法分为基于灰度的和基于特征的配准方法。

基于灰度的方法通常是通过优化3D术前图像和2D术中图像投影的相似性测度来实现的。数字重建放射成像(DRR)和最大强度投影(MIP)是两种常见的通过计算机断层成像(CT)图像生成模拟X射线投影的方法。Hipwell分析了DRR/MIP与X射线图像的六种相似度指标，其中模式强度和梯度差异表现最好。考虑到基于灰度的方法使用整个图像灰度信息进行配准，因此它们对背景异常值非常敏感。此外，当配准具有大尺度变换的数据时，基于优化的方法可能产生较小的捕获范围。

基于特征的3D/2D配准依赖于从两种模态的图像中提取的具有一致性特征。中心线是血管配准中最常用的特征表示。迭代最近点(ICP)方法将点云配准分解为包含匹配和配准阶段的交替连续的过程。对于匹配阶段，可以通过查找点的最小欧氏距离来分配点对应。Baka介绍了一种将二维点反投影到三维空间，然后执行经典ICP过程的3D/2D配准方法。Rivest-Henault利用血管中心预先计算距离变换，并构造目标函数，可以加速配准过程。Benseghir提出了迭代最近邻曲线(ICC)方法，该方法利用血管分支作为基于最近邻关系的配对元素，然后估计出最小化配对分支之间距离和的变换。这些类ICP方法对噪声和异常值很敏感，因为它们将对应关系限制为一对一的匹配。基于最近邻关系的匹配也会导致这些方法严重依赖于初始姿态。

由于硬分配策略对噪声和异常值的敏感性，在概率分配框架中，软分配策略将一对一对应关系松弛为一对多。在高斯混合模型(GMM)和期望最大化(EM)算法的基础上，Myronenko提出了相干点漂移(CPD)方法，该方法强制GMM中心体进行相干移动，以保持点集的全局拓扑结构。Kang使用与CPD相同的框架进行3D/2D点云配准。在考虑透视投影的非线性特性时，采用粒子群优化算法求解配准参数的最优估计。Baka提出了一种从Jian andVemuri的方法扩展而来的方向约束OGMM方法，利用方向和位置估计两个点集的L₂距离，并对其进行优化，实现3D/2D冠状动脉配准。由于OGMM利用了中心线的定位，因此对噪声数据具有较高的精度和鲁棒性。

由于血管拓扑是跨形态和维度的不变属性，因此图匹配成为一种有效的血管配准方法。血管图的匹配可以描述为先估计血管分叉点的对应关系，然后将分叉点视为顶点来匹配连接它们之间的曲线。Serradell将血管配准描述为寻找最大可能对应关系的搜索过程，并使用先验搜索来加速该过程。

Pinheiro基于待匹配图的拓扑一致性，将血管匹配描述为一种树搜索方法，并采用蒙特卡罗树搜索来解决该问题。Moriconi通过设置节点和边属性来定义图匹配的亲密度函数，然后将二次分配问题的函数最大化，得到节点匹配。

上述几个方法均与3D/3D或2D/2D血管系统匹配/配准相关。对于3D/2D血管匹配，投影中的重叠问题决定了该方法需要强调异常值和噪声。ICC方法通过两个节点之间的配对曲线来使用血管拓扑。对于有噪声的二维图，利用邻域关系可以约束候选对象。通过保证分叉点的连通性，Benseghir在ICC方法的基础上提出了一种分而治之的树形中心线匹配方法。Liu将3D和2D血管视为树形拓扑，将树形拓扑表示为序列。然后，通过序列对序列的匹配，实现节点的血管匹配。通过保持前后关系，利用拓扑排序算法提取序列，然后按顺序遍历。

由于透视投影的非线性性质，基于数值优化的3D/2D配准容易陷入局部极值，使得这些方法对初始配准位姿敏感。初始化操作是3D/2D配准的关键。一般来说，术前和术中图像来自不同的设备。大多数配准方法的捕获范围不足以覆盖两种采集设备的坐标系统之间的转换。Markelj介绍了一些初始化方法，如患者位置和方向的对齐、相应标记对的注册以及手动初始化。然而，一种利用内在特征的自动初始化方法可能更适用于术中3D/2D的配准。Varnavas预先计算了大范围3D位姿的二维投影模板，并使用术中图像评估模板与二维透视的广义霍夫变换之间的相似性，以获得初始对准。Miao建立了从金属植入体轮廓中提取的二维轮廓的形状上下文编码库，并采用Jensen-Shannon divergence算法作为快速库匹配的匹配度量。Gouveia提出了一种基于回归的CTA和XRA的初始配准方法，该方法将二维投影图像特征与刚性变换参数联系起来。

为了获得较大的捕获范围或对初始位姿不敏感的配准方法，基于“先匹配再求变换”的框架比基于优化的方法更合适。假定已知3D和2D点的对应关系，3D/2D配准的目标类似于计算机视觉领域的PnP(perspective-n-point)问题。它的目的是通过将对点变换的估计问题简化为对四个控制点坐标的估计问题来确定摄像机的姿态。

发明内容

为克服现有技术的缺陷，本发明要解决的技术问题是提供了一种基于启发式树搜索的血管结构3D/2D刚性配准方法，其3D图像和2D图像的血管图匹配精确度高。

本发明的技术方案是：这种基于启发式树搜索的血管结构3D/2D刚性配准方法，其包括以下步骤：

(1)利用3D和2D血管拓扑一致性来实现血管图匹配，将图的匹配结果表示为成对的3D和2D边的集合；

(2)将血管匹配过程看作是在已有的血管对上增加一对新的匹配边的连续过程，将匹配过程分解为连续状态，并用于搜索树的构造；

(3)在搜索树的每个节点上，使用一个封闭解来计算基于血管点密集匹配的配准结果，设计一个评价配准质量和配准质量的节点评分；

(4)在A-star搜索算法的基础上，利用一种改进的启发式树搜索策略，寻找节点得分最高的最优结果。

本发明利用3D和2D血管拓扑一致性来实现血管图匹配，图的匹配结果可以表示为成对的3D和2D边的集合，根据血管中心线的拓扑连续性，血管匹配过程可以看作是在已有的血管对上增加一对新的匹配边的连续过程，此属性可将匹配过程分解为连续状态，并用于搜索树的构造，将血管图的3D/2D配准问题表示为一个树搜索问题，在树的每个节点上，使用一个封闭解来计算基于血管点密集匹配的配准结果，设计了一个评价配准质量和配准质量的节点评分，在A-star搜索算法的基础上，提出了一种改进的启发式树搜索策略，寻找节点得分最高的最优结果，从而使得3D图像和2D图像的血管图匹配精确度高。

还提供了基于启发式树搜索的血管结构3D/2D刚性配准装置，其包括：

匹配结果表示模块，其利用3D和2D血管拓扑一致性来实现血管图匹配，将图的匹配结果表示为成对的3D和2D边的集合；

搜索树的构造模块，其将血管匹配过程看作是在已有的血管对上增加一对新的匹配边的连续过程，将匹配过程分解为连续状态，并用于搜索树的构造；

计算结果及节点评分模块，其在搜索树的每个节点上，使用一个封闭解来计算基于血管点密集匹配的配准结果，设计一个评价配准质量和配准质量的节点评分；

最优结果寻找模块，其在A-star搜索算法的基础上，利用一种改进的启发式树搜索策略，寻找节点得分最高的最优结果。

附图说明

图1a是待配准的3D血管模型；图1b是3D血管中心线；图1c是2D造影图像；图1d是2D血管中心线。

图2是根据本发明的血管匹配的部分搜索树。

图3示出了根据本发明的基于启发式树搜索的血管结构3D/2D刚性配准方法的流程图。

具体实施方式

如图3所示，这种基于启发式树搜索的血管结构3D/2D刚性配准方法，其包括以下步骤：

(1)利用3D和2D血管拓扑一致性来实现血管图匹配，将图的匹配结果表示为成对的3D和2D边的集合；

(2)将血管匹配过程看作是在已有的血管对上增加一对新的匹配边的连续过程，将匹配过程分解为连续状态，并用于搜索树的构造；

(3)在搜索树的每个节点上，使用一个封闭解来计算基于血管点密集匹配的配准结果，设计一个评价配准质量和配准质量的节点评分；

(4)在A-star搜索算法的基础上，利用一种改进的启发式树搜索策略，寻找节点得分最高的最优结果。

本发明利用3D和2D血管拓扑一致性来实现血管图匹配，图的匹配结果可以表示为成对的3D和2D边的集合，根据血管中心线的拓扑连续性，血管匹配过程可以看作是在已有的血管对上增加一对新的匹配边的连续过程，此属性可将匹配过程分解为连续状态，并用于搜索树的构造，将血管图的3D/2D配准问题表示为一个树搜索问题，在树的每个节点上，使用一个封闭解来计算基于血管点密集匹配的配准结果，设计了一个评价配准质量和配准质量的节点评分，在A-star搜索算法的基础上，提出了一种改进的启发式树搜索策略，寻找节点得分最高的最优结果，从而使得3D图像和2D图像的血管图匹配精确度高。

优选地，所述步骤(1)中，血管图表示为其中/>为顶点集合，ε为边的集合，/>包含血管的端点和分叉点，边表示为连接两个相邻顶点的曲线，3D血管图/>和2D血管图/>的匹配分别表示为顶点对的集合/>边对的集合/>π^ε通过/>直接导出。

优选地，所述步骤(1)中，根据公式(1),3D和2D血管结构的配准表示为找到3D血管的最优变换T，在该变换T作用后的3D血管的投影与2D血管达到最佳对齐，

其中T表示3D刚性变换，表示两个血管图的量化距离，P为透视投影操作，它是由X射线成像设备决定的且不变的；根据公式(2),当已知血管边的匹配时，配准重新描述为找到3D血管的最优变换，使每个3D血管边与2D血管边达到最佳对齐

因此，配准问题则转变成同时计算变换T和匹配π^ε的问题。

优选地，所述步骤(2)中，根据π^ε的表达式，血管图的匹配描述为逐步增加新的匹配边对的连续过程这个特性将匹配过程分解为连续的状态，并且利用这些状态构建搜索树；在构建搜索树的过程中，为了保证拓扑连续性，新加入的血管边对/>满足以下两个条件：

(a)和/>的开始顶点包含于已匹配的顶点集合中/>而它们的结束顶点不在其中；

(b)不存在两条边重叠的情况，这个步骤表达为有效配对探测，给定两个血管图和/>以及当前状态的匹配结果/>获得多组候选匹配边/>

优选地，图2中展示了一个搜索树的部分示例，在树的根节点选中一对顶点，并沿着树开始拓展分支。所述步骤(2)中，在一个搜索树的根节点选中一对顶点，并沿着树开始拓展分支，每一次拓展的新的树节点，边的匹配结果πε更新；节点的匹配结果由π^ε决定，也更新；新的刚性变换T通过解决公式(2)来获得；一个于匹配π^ε和配准T相关的分数/>被计算得到，搜索树的每个节点定义为/>

优选地，所述步骤(3)中，

假设边的匹配结果π^ε给定，刚性变换结果通过闭合解T求得；对于3D点云C^3D和2D点云C^2D，配准问题重新定义为公式(3)

使用闭合解的方式来获得刚性变换，将3D点的坐标用在世界坐标系下的四个控制点/>的加权来表示，/>满足/>其中世界坐标系用上标/>表示，相机坐标系用上标^(c)表示，/>的坐标信息通过系数α_ij和虚拟控制点获得，求取由世界坐标系到相机坐标系中3D点的转换映射T:/>等同于虚拟控制点对的映射T:/>因此，通过计算相机坐标系中控制点{u_j ^(c)}坐标，获取刚性变换矩阵T；不同坐标系下，虚拟控制点和3D血管点具有相同的系数α_ij，因此，利用给定的矩阵π，假设每个3D点投影后的坐标与其对应的2D点的坐标一致，得到/>由于/>对上式进行重写，得到公式(4)

其中H是相机内部校准矩阵，ω_i是标量投影参数,通过求解公式

(4)得到相机坐标系下的四个虚拟控制点的坐标，刚性变换矩阵T:则利用对应点对映射关系计算得到。

优选地，所述步骤(3)中，为了评估搜索树中每个节点的配准和匹配结果，定义一个与π^ε和T有关的分数

优选地，所述步骤(3)中，包含两个准则：

第一个准则对应于一个期望，鼓励更多的具有高精确度的成对边，

第一个准则的分数定义为公式(5)

其中d_F(·)表示两个点集序列的Fréchet距离，σ是用于归一化该距离的尺度参数，与π^ε和T相关，它的上确界是匹配边对的数量|π^ε|；第二条准则对应一个通用的配准假设，越多3D点的投影与2D点重合，则配准结果越好，/>通过公式(6)获得

其中|C^3D|表示点集C^3D的数量，的上确界是2，3D点的投影/>

和2D点/>的最短距离通过距离变换计算得到。

优选地，所述步骤(3)中，

总体的分数表示为公式(7)

其中α是用于平衡和/>贡献的系数。

HTSR方法的伪代码如下所示。

本领域普通技术人员可以理解，实现上述实施例方法中的全部或部分步骤是可以通过程序来指令相关的硬件来完成，所述的程序可以存储于一计算机可读取存储介质中，该程序在执行时，包括上述实施例方法的各步骤，而所述的存储介质可以是：ROM/RAM、磁碟、光盘、存储卡等。因此，与本发明的方法相对应的，本发明还同时包括一种基于启发式树搜索的血管结构3D/2D刚性配准装置，该装置通常以与方法各步骤相对应的功能模块的形式表示。该装置包括：

匹配结果表示模块，其利用3D和2D血管拓扑一致性来实现血管图匹配，将图的匹配结果表示为成对的3D和2D边的集合；

搜索树的构造模块，其将血管匹配过程看作是在已有的血管对上增加一对新的匹配边的连续过程，将匹配过程分解为连续状态，并用于搜索树的构造；

计算结果及节点评分模块，其在搜索树的每个节点上，使用一个封闭解来计算基于血管点密集匹配的配准结果，设计一个评价配准质量和配准质量的节点评分；

最优结果寻找模块，其在A-star搜索算法的基础上，利用一种改进的启发式树搜索策略，寻找节点得分最高的最优结果。

以上所述，仅是本发明的较佳实施例，并非对本发明作任何形式上的限制，凡是依据本发明的技术实质对以上实施例所作的任何简单修改、等同变化与修饰，均仍属本发明技术方案的保护范围。

Claims (6)

1.基于启发式树搜索的血管结构3D/2D刚性配准方法，其特征在于：其包括以下步骤：

(1)利用3D和2D血管拓扑一致性来实现血管图匹配，将图的匹配结果表示为成对的3D和2D边的集合；

(2)将血管匹配过程看作是在已有的血管对上增加一对新的匹配边的连续过程，将匹配过程分解为连续状态，并用于搜索树的构造；

(3)在搜索树的每个节点上，使用一个封闭解来计算基于血管点密集匹配的配准结果，设计一个评价配准质量和配准质量的节点评分；

(4)在A-star搜索算法的基础上，利用一种改进的启发式树搜索策略，寻找节点得分最高的最优结果；

所述步骤(3)中，

假设边的匹配结果π^ε给定，刚性变换结果通过闭合解T求得；对于3D点云C^3D和2D点云C^2D，配准问题重新定义为公式(3)

使用闭合解的方式来获得刚性变换，将3D点的坐标用在世界坐标系下的四个控制点/>的加权来表示，/>满足/>其中世界坐标系用上标/>表示，相机坐标系用上标^(c)表示，/>的坐标信息通过系数α_ij和虚拟控制点获得，求取由世界坐标系到相机坐标系中3D点的转换映射T:/>等同于虚拟控制点对的映射T:/>因此，通过计算相机坐标系中控制点{u_j ^(c)}坐标，获取刚性变换矩阵T；不同坐标系下，虚拟控制点和3D血管点具有相同的系数α_ij，因此，利用给定的矩阵π，假设每个3D点投影后的坐标与其对应的2D点的坐标一致，得到/>由于/>对上式进行重写，得到公式(4)

其中H是相机内部校准矩阵，ω_i是标量投影参数,通过求解公式(4)得到相机坐标系下的四个虚拟控制点的坐标，刚性变换矩阵T:则利用对应点对映射关系计算得到；

所述步骤(3)中，为了评估搜索树中每个节点的配准和匹配结果，定义一个与π^ε和T有关的分数

所述步骤(3)中，包含两个准则：

第一个准则对应于一个期望，鼓励更多的具有高精确度的成对边，第一个准则的分数定义为公式(5)

其中d_F(·)表示两个点集序列的Fréchet距离，σ是用于归一化该距离的尺度参数，与π^ε和T相关，它的上确界是匹配边对的数量|π^ε|；第二条准则对应一个通用的配准假设，越多3D点的投影与2D点重合，则配准结果越好，/>通过公式(6)获得

其中|C^3D|表示点集C^3D的数量，的上确界是2，3D点的投影/> 和2D点/>的最短距离通过距离变换计算得到；

所述步骤(3)中，

总体的分数表示为公式(7)

其中α是用于平衡和/>贡献的系数。

2.根据权利要求1所述的基于启发式树搜索的血管结构3D/2D刚性配准方法，其特征在于：所述步骤(1)中，血管图表示为 其中/>为顶点集合，ε为边的集合，/>包含血管的端点和分叉点，边表示为连接两个相邻顶点的曲线，3D血管图/>和2D血管图/>的匹配分别表示为顶点对的集合/>边对的集合/>π^ε通过/>直接导出。

3.根据权利要求2所述的基于启发式树搜索的血管结构3D/2D刚性配准方法，其特征在于：所述步骤(1)中，根据公式(1),3D和2D血管结构的配准表示为找到3D血管的最优变换T，在该变换T作用后的3D血管的投影与2D血管达到最佳对齐，

其中T表示3D刚性变换，表示两个血管图的量化距离，P为透视投影操作，它是由X射线成像设备决定的且不变的；根据公式(2),当已知血管边的匹配时，配准重新描述为找到3D血管的最优变换，使每个3D血管边与2D血管边达到最佳对齐

因此，配准问题则转变成同时计算变换T和匹配π^ε的问题。

4.根据权利要求3所述的基于启发式树搜索的血管结构3D/2D刚性配准方法，其特征在于：所述步骤(2)中，根据π^ε的表达式，血管图的匹配描述为逐步增加新的匹配边对的连续过程 这个特性将匹配过程分解为连续的状态，并且利用这些状态构建搜索树；在构建搜索树的过程中，为了保证拓扑连续性，新加入的血管边对满足以下两个条件：

(a)和/>的开始顶点包含于已匹配的顶点集合/>中，而它们的结束顶点不在其中；

(b)不存在两条边重叠的情况，这个步骤表达为有效配对探测，给定两个血管图和以及当前状态的匹配结果/>获得多组候选匹配边/>

5.根据权利要求4所述的基于启发式树搜索的血管结构3D/2D刚性配准方法，其特征在于：所述步骤(2)中，在一个搜索树的根节点选中一对顶点，并沿着树开始拓展分支，每一次拓展的新的树节点，边的匹配结果π^ε更新；节点的匹配结果由π^ε决定，也更新；新的刚性变换T通过解决公式(2)来获得；一个于匹配π^ε和配准T相关的分数/>被计算得到，搜索树的每个节点定义为/>

6.基于启发式树搜索的血管结构3D/2D刚性配准装置，其特征在于：其包括：

匹配结果表示模块，其利用3D和2D血管拓扑一致性来实现血管图匹配，将图的匹配结果表示为成对的3D和2D边的集合；

搜索树的构造模块，其将血管匹配过程看作是在已有的血管对上增加一对新的匹配边的连续过程，将匹配过程分解为连续状态，并用于搜索树的构造；

计算结果及节点评分模块，其在搜索树的每个节点上，使用一个封闭解来计算基于血管点密集匹配的配准结果，设计一个评价配准质量和配准质量的节点评分；

最优结果寻找模块，其在A-star搜索算法的基础上，利用一种改进的启发式树搜索策略，寻找节点得分最高的最优结果；

所述计算结果及节点评分模块中，

假设边的匹配结果π^ε给定，刚性变换结果通过闭合解T求得；对于3D点云C^3D和2D点云C^2D，配准问题重新定义为公式(3)

使用闭合解的方式来获得刚性变换，将3D点的坐标用在世界坐标系下的四个控制点/>的加权来表示，/>满足/>其中世界坐标系用上标/>表示，相机坐标系用上标(c)表示，/>的坐标信息通过系数α_ij和虚拟控制点获得，求取由世界坐标系到相机坐标系中3D点的转换映射T:/>等同于虚拟控制点对的映射T:/>因此，通过计算相机坐标系中控制点{u_j ^(c)}坐标，获取刚性变换矩阵T；不同坐标系下，虚拟控制点和3D血管点具有相同的系数α_ij，因此，利用给定的矩阵π，假设每个3D点投影后的坐标与其对应的2D点的坐标一致，得到/>由于/>对上式进行重写，得到公式(4)

其中H是相机内部校准矩阵，ω_i是标量投影参数,通过求解公式(4)得到相机坐标系下的四个虚拟控制点的坐标，刚性变换矩阵T:则利用对应点对映射关系计算得到；

所述步骤(3)中，为了评估搜索树中每个节点的配准和匹配结果，定义一个与π^ε和T有关的分数

所述计算结果及节点评分模块中，包含两个准则：

第一个准则对应于一个期望，鼓励更多的具有高精确度的成对边，第一个准则的分数定义为公式(5)

其中d_F(·)表示两个点集序列的Fréchet距离，σ是用于归一化该距离的尺度参数，与π^ε和T相关，它的上确界是匹配边对的数量|π^ε|；第二条准则对应一个通用的配准假设，越多3D点的投影与2D点重合，则配准结果越好，/>通过公式(6)获得

其中|C^3D|表示点集C^3D的数量，的上确界是2，3D点的投影/> 和2D点/>的最短距离通过距离变换计算得到；

所述计算结果及节点评分模块中，

总体的分数表示为公式(7)

其中α是用于平衡和/>贡献的系数。