CN101763642A - 一种用于冠脉造影三维重建的匹配方法 - Google Patents

Abstract

一种用于冠脉造影三维重建的匹配方法，属于数字图像处理与医学成像的交叉领域，目的是满足临床医学上心血管疾病的辅助检测以及手术导航的特殊要求。本发明包括模型指导造影图血管段匹配以及造影图血管点匹配两个步骤。本发明使用多尺度的血管树模型，更好地指导血管匹配，同时可以很好的解决个体差异性的问题，通过迭代匹配及搜索势能最大下降方法，提高了血管段和血管点的匹配精度。本发明可以得到很好的血管造影匹配结果，解决了由多视角血管造影图进行可靠的自动三维重建难题。

Description

一种用于冠脉造影三维重建的匹配方法

技术领域

本发明属于数字图像处理与医学成像的交叉领域，具体涉及一种模型指导下的血管造影匹配方法。

背景技术

血管树三维重建是通过不同视角的X射线二维投影图像中相应的图像信息恢复血管三维空间结构的过程。它一般包括3个步骤：(1)血管骨架的提取；(2)血管段的识别和匹配；(3)血管空间点的估计和重建及血管基元的拟合。而在这三个步骤中，最关键同时也是最难的一步就是两幅造影图中各血管段的识别及其上各血管点的匹配。

对于心血管双臂造影图像的结构识别和匹配问题，存在以下几个难点：第一，两个造影成像面中表现的血管拓扑结构往往相差较大；第二，投影图像中血管投影交叉的存在，严重影响着血管拓扑结构的识别和描述；第三，造影过程中血管段经常会发生遮挡使得大量信息发生缺失，不利于血管段上各血管点的匹配。

传统的方法是用通过外极线约束来找两幅造影图中匹配的血管点对，但是外极线约束极易发生误匹配，需要有其他一些约束条件或知识来指导造影图中血管段结构的识别与匹配。现有技术中关于基于双臂造影系统的常见方法有：

1、基于知识指导的匹配方法。M.Garreau等人提出了利用结构学、拓扑学以及解剖学等知识对每幅视图中血管进行标记，来指导血管的三维重建。G.Coppini等人为实现三维血管骨架的重建，使用了血管的几何特征、拓扑结构以及解剖学等知识来指导两幅造影图中血管骨架的匹配，然后将重建得到的血管骨架作反投影，计算其误差，校验并优化匹配结果。D.Delaere等人提出了一种基于知识的三维重建系统，经过外极线约束，相异性约束以及解剖学约束(即采用左冠状动脉(LCA)解剖学静态模型)，对两幅在标准投影角度下获得的造影图进行匹配。

基于知识指导的血管段匹配方法实现比较简单，但是由于知识指导太过抽象，它很难对各个体都取得很好的效果，同时由于其中用到了血管的结构和几何信息，对造影图的质量也有较高要求。

2、基于拓扑相似性的匹配方法：I.Liu等人设计了一种对互成90°角的两幅血管投影图进行三维重建的算法，通过外极线约束，血管段的几何相似性和连通性来实现两幅造影图中血管段的匹配。黄家祥提出综合利用拓扑约束和外极线约束匹配二维骨架树的分支点和端点，从而实现两幅造影图上血管段的匹配，在造影图像质量较好时能达到较高的匹配精度。但是当造影质量较差，从两幅造影图中分割得到的血管拓扑结构有较大区别时匹配效果就不是很好。

基于拓扑相似性的血管段匹配方法对于两幅造影图的成像质量要求都很高，在成像效果好，即两幅造影图中提取得到的血管拓扑结构相同时，能得到较好的匹配效果；但是当两幅图中有一幅出现血管间的遮挡、重叠或血管段的缺失时，则很难得到正确的匹配关系。

3、基于蛇形的匹配方法

Radeva等人提出了一种基于能量最小化的三维重建算法。将三维血管用蛇形表示，利用血管的平滑性、连通性以及空间结构(内力)和三维重建结果的反投影误差(外力)构建能量函数。通过蛇形变化获得两图中血管的映射关系，并通过优化系统中心(isocenter)的位置来减少匹配误差。Cristina

等人结合预测的(非)扭曲模型((un)distortion model)，利用双平面蛇形技术(BiplaneSnake)对血管造影图中的血管中轴线进行三维重建，并应用一般梯度矢量流(GGVF)优化蛇形波浪式运动，使重建结果更加完善。

基于蛇形的血管匹配方法需要进行人工干预，要预先选好初始匹配点对，而且构造能量函数的过程太过复杂，实现起来较难。

4、基于个体模型数据库的匹配方法：

P.Hall等人基于个体模型数据库的匹配方法提出了一种具有学习功能的系统，通过存储包括各种变化和个体差异性的血管数据建立一个血管模型数据库。首先根据学习系统中血管目录(VC)与数据造影图的映射关系，产生假设，然后参照各对应血管段的具体特征进行假设验证，继而获得重建结果。基于个体模型数据库利用的是个体冠状动脉之间的相似性进行匹配，如果模型库够大可能得到较好的匹配效果，但是基于个体模型库的方法存在一个数据量和存储的问题。要得到这么大量的数据比较困难，即使能得到也需要有一个很大的存储空间来存储这些数据，工程实在太过庞大，且具体算法实现时，计算复杂，耗时长。

5、基于血管模型的匹配方法：

C.Smets等人介绍了一种基于血管解剖学模型的匹配方法，该模型是基于规则的，它由一些标准参考资料和心脏病专家提出的约束条件组成。该方法首先采用一般血管模型描述各血管段，然后通过解剖学模型中的约束条件进行标记，最后进行匹配和三维重建。Claire Chalopin等人提出了一种用血管树模型对造影图中血管进行标记的方法，利用Renaudin等人建立的一个通用冠状动脉模型Coronix，将模型的投影图与造影图进行匹配从而分别对造影图中的血管按模型中的血管名字进行命名，然后可以通过命名结果对不同角度的造影图中的血管进行匹配。利用静态通用模型指导造影图中血管标记的方法克服了基于知识指导过于抽象的缺点，算法直观明白，能得到一个比较稳定的效果。但是它过多依靠模型的具体数值，对于个体特异性的问题无法得到很好的解决。

在中国专利文献CN101301207A中也公开了一种用于冠脉造影三维重建的匹配方法，但是，其在进行血管段匹配时，仅利用了模型投影图血管主子树与造影图血管子树的拓扑结构相似程度作为匹配评价，且对血管点的匹配，也只是在匹配血管段上等间隔取样来选取匹配点进行三维重建，这样做重建误差会较大，因为在出现血管缺失、重叠等情况下，并不能保证左右造影图中的血管段和取样点是一一对应的关系。

本发明提出的用于冠脉造影三维重建的匹配方法，是对上述专利的进一步深化与改进，并相对于之前所介绍的常见匹配方法，在图像要求质量高，耗时长，需人工干预多及由个体差异引起的引起鲁棒性差等方面都有了明显的改善。实验证明，该方法不仅对双臂造影图像可以得到很好的匹配效果，对于单臂造影图像的匹配也能得到较高的匹配精度。该方法可以解决由多视角血管造影图进行可靠的自动三维重建难题，满足临床医学心血管疾病辅助检测和手术导航的应用要求。

发明内容

本发明提供一种模型指导的血管造影匹配方法，通过心脏与心血管的静态模型，利用模型投影图血管主子树与造影图血管子树的拓扑结构相似程度以及两幅不同视角的X射线血管造影图(即左右造影图)的重建反投影误差作为匹配评价，指导双视角造影图中的血管段匹配，然后通过重建反投影误差构成的势能梯度下降方法获得局部最优的匹配点来实现血管点匹配。

一种用于冠脉造影三维重建中的匹配方法，包括：

(一)、造影图血管段的匹配，即对从左右造影图中提取的血管中轴图像中的血管段进行匹配，主要包括：

(1)提取基础数据：从两幅不同视角的X射线血管造影图中，分别提取血管中轴图像，从而获得各自的血管树即血管拓扑结构，及血管的半径信息；

(2)预匹配：通过血管半径信息识别并匹配左右造影图中左冠脉或右冠脉的两个主要分支，确定所述两个主要分支各自的主子树，从而确定所述主子树在左右造影图中的对应关系；

(3)对所述主子树进行迭代识别和匹配，具体包括以下两个步骤：

A.模型指导下的部分匹配步骤，即按X射线造影方式从不同角度对血管树静态模型进行投影，通过投影图与造影图之间的匹配对左右两幅造影图中部分血管进行分别命名标记，从而实现左右造影图中部分血管段的匹配，具体过程为：搜索左右造影图的两个主子树与模型投影图的相对应的主子树之间所有对应关系，利用最小化匹配代价原则找到最佳投影角度与最优匹配方式，即在造影图的血管树上找到与模型投影图拓扑结构最相似的主子树，进而对造影图中部分血管段进行标记，实现部分匹配；

B.主子树其余部分的匹配步骤，根据拓扑结构的约束和最小化反投影误差对它们进行迭代匹配，具体匹配过程如下：

设左右两幅造影图中血管拓扑结构分别表示为left_T和right_T，对于所述两幅造影图中主子树上还没有命名标记的所述其余部分血管段，存在三种对应情况：

第一种情况，left_T中未命名血管段UnName1在right_T中没有对应血管段，则对UnName1不予处理；

第二种情况，left_T中未命名血管段UnName1在right_T中恰好有一个对应血管段UnName2，则令UnName1和UnName2为匹配血管，若UnName1和UnName2不止是一段血管，下面还有分支的话，则通过搜索left_T与right_T中UnName1和UnName2及它们的分支之间的所有对应关系，利用最小化反投影误差得到最优匹配方式；

第三种情况，left_T中未命名血管段UnName1在right_T中有两个以上对应的血管段，那么未命名血管段UnName1有可能与其中任何一个血管段匹配，因此需分别计算所有可能情况下的代价函数g，取其中代价函数最小的为匹配血管段，代价函数g通过下式计算：

g＝Err_pro(I1_1)+Err_pro(I1_2)+Err_pro(UnName1，UnName2)

其中，UnName2为所述两个以上对应的血管段中任一个血管段，I1_1和I1_2表示left_T(或right_T)中标记为I1_1和I1_2的两个首尾相连的血管段，在left_T中，I1_1和I1_2被UnName1所分，在right_T中，I1_1和I1_2被UnName2所分，Err_pro(I1_1)和Err_pro(I1_2)表示分别将左右造影图中标记为I1_1和I1_2的血管段进行重建的反投影误差，Err_pro(UnName1，UnName2)表示将UnName1和UnName2进行重建的反投影误差；

(二)、造影图血管点匹配，即对匹配好的血管段中的各血管点进行匹配，包括以下步骤：

(1)首先通过外极线约束得到一系列初始的血管匹配点对(P′_i，Q′_i)(i＝1，2，...，m，m为血管匹配点对的个数)；

(2)由于匹配点在血管段上的空间位置需相对应，基于此顺序约束，调整初始匹配点对，得到调整后的血管匹配点对(P_i，Q_i)；

(3)找到调整后的血管匹配点对(P_i，Q_i)中势能梯度最大的点对，在保持所述顺序约束条件下，沿着血管以像素为单位向势能降低最大的方向调整此匹配点对的位置；

其中，匹配点对(P_i，Q_i)的势能即为对这两个点进行三维重建后的反投影误差Err(P_i，Q_i)，而匹配点对(P_i，Q_i)的势能梯度定义为：

Grad(P_i，Q_i)＝max((Err(P_i，Q_i)-Err(P_i，Q_i_pre))，(Err(P_i，Q_i)-Err(P_i，Q_i_next))，

(Err(P_i，Q_i)-Err(P_i_pre，Q_i))，(Err(P_i，Q_i)-Err(P_i_next，Q_i)))

上式中，P_i_pre，P_i_next，Q_i_pre和Q_i_next分别为点P_i和Q_i邻域的前点与后点，且它们均在血管段上；

(4)重复上述步骤(3)，直到所有匹配点对都达到一个最低势能状态。

进一步地，最小化匹配代价原则具体为，对于每个投影角度下所有满足对应关系的情况，比较其中每种匹配方式的匹配代价f₁和f₂，f₁为左右造影图之间的匹配代价，f₂为模型投影图与造影图之间的匹配代价。当f₁具有唯一最小值时，则该f₁所对应的投影角度和匹配方式为最佳投影角度与最优匹配方式，当f₁最小值不唯一时，那么则继续比较在这些情况下的f₂，取f₂最小时对应的投影角度和匹配方式为最佳投影角度与最优匹配方式，然后通过最佳投影角度与最优匹配方式利用模型对造影图中部分血管进行标记。

其中，匹配代价f₁和f₂的具体计算方法如下：

(1)左右造影图主子树之间的匹配代价f₁的计算

取f₁为对左右造影图中匹配血管段进行重建后的反投影与原血管段的误差的平均值：

$f_1=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\sqrt{\mathrm{dis}{(P_{1i}-Q_{1i})}^2+\mathrm{dis}{(P_{2i}-Q_{2i})}^2}/N$

其中P_1i和P_2i分别为对应血管段上选取的对应血管点，Q_1i和Q_2i是将P_1i和P_2i进行重建之后按左右造影图参数反投影的点，N是选取的点对的个数。

(2)模型投影图主子树与造影图主子树之间的匹配代价f₂的计算

f₂＝ω₁*f(angle)+ω₂*f(length)

$f\left(\mathrm{angle}\right)=\underset{i=1}{\overset{M}{\mathrm{\Σ}}}|{\mathrm{angle}}_{1i}-{\mathrm{angle}}_{2i}|$

$f\left(\mathrm{length}\right)=\underset{j=1}{\overset{L}{\mathrm{\Σ}}}\sqrt{|l_{11j}*l_{22j}-l_{12j}*l_{21j}|}$

其中f(angle)表示两者之间血管拓扑结构的匹配代价，angle_1i和angle_2i分别为模型投影图和造影图中对应血管段的第一个点与其中心点连接所成的角度，f(length)表示两段血管长度的匹配代价，l_11j和l_12j，l_22j和l_21j分别为模型投影图和造影图中对应分叉点的左右子血管段的长度，ω₁和ω₂分别为f(angle)和f(length)的匹配代价系数，这里令ω₁＝mean(l_2j)，ω₂＝1，M是模型投影图与造影图对应血管段的数目，L是模型投影图与造影图对应分叉点的数目。

进一步地，搜索对应关系的具体过程为：

给造影图中血管树上任意一个分叉点设计一个信号量S_k，约定当S_k＝0时，该分叉点在子树中也是一个分叉点；当S_k＝1时，该分叉点的父血管段与左子血管段在子树中连接成一段血管；当S_k＝2时，该分叉点的父血管段与右子血管段在子树中连接成一段血管。设造影图中血管树上有n个分叉点，对S_k(k＝1，2，...，n)进行3进制编码，则共存在3ⁿ种编码情况，每一种编码情况对应一种可能与模型投影图中血管树拓扑结构相似的血管结构。对每一种编码，判断其对应的血管结构是否与模型投影图中血管树的拓扑结构相似，如果相似即确定了两者是对应的。

对每种对应情况，计算左右造影图之间的匹配代价f₁，以及模型投影图主子树与造影图主子树之间的匹配代价f₂。

本发明方法考虑到人体冠状动脉系统结构的拓扑相似性，用心血管树模型分别对两个不同视角的X射线造影图中的血管基元投影进行识别，从而指导两幅造影图中血管基元投影的匹配，提高血管结构特征匹配的可靠性和准确度，增大血管三维重建的精度。具体而言，本发明具有如下四个方面的技术效果：

(1)使用多尺度的血管树模型，更好地指导血管匹配。国外也有用模型指导血管标记与匹配的算法，但他们使用的模型都是一个固定的模型，没有尺度级的变化，不能满足不同个体对于细节的要求。这里我们在通过用模型指导双视角造影图匹配的过程中，同时对模型进行补充，可添加以前没有的细节血管。

我们提出与血管树级数等价的多尺度的血管树模型的概念，血管树多尺度是相对血管树的级数来说的。根据血管树的解剖学知识，血管树可以分成很多不同级数的血管，第一级血管是最根部的血管，在冠脉系统中就是连接主动脉的冠状动脉，它的分支就是二级血管，同理，三级血管是二级血管的分支，依次往下。在重建的时候，我们可以根据重建要求和造影图的实际情况，选择使用不同尺度或不同级数的血管树模型对造影图进行匹配指导，从而可以更有效地指导血管的配准与重建。比如说重建只要求到第三级血管，我们就只需要提供三级血管树模型指导匹配。或者分割后造影图中能得到的最大细节只能达到三级，我们同样只需要用三级血管树模型来指导匹配。通过建立多尺度的血管树模型提高了血管匹配和三维重建的可控性和鲁棒性，提高了重建精度。

(2)血管匹配算法并没有过分依赖于模型几何信息的具体数值，可以很好的解决个体差异性的问题。在国外的很多算法中，用模型指导血管匹配都利用了很多关于模型几何的信息，比如投影图中血管的坐标，血管基元长度等，而这些会受到血管尺寸、扭曲、以及偏移等个体差异性的影响，把它用在算法中会大大降低算法的鲁棒性；相比之下，拓扑结构就比几何信息稳定得多，它不会受到图像尺寸的影响，同时具有旋转不变性，且不受平移影响。鉴于此，我们的方法只利用了血管的拓扑结构和部分对图像尺寸和旋转不太敏感的几何信息，如相邻血管基元长度比，血管基元曲率(定义为在血管基元上所有点的平均曲率)等，来指导血管的匹配，从而使本方法能降低对个体差异性的敏感性，使重建过程变得更加鲁棒和更高自动化程度。

(3)在利用模型拓扑结构对造影图进行匹配的基础上，通过迭代匹配获取重建的反投影误差最小的匹配结果作为最终的血管段匹配结果，大大提高了血管段的匹配精度，既利用了冠脉模型的拓扑结构特征，同时很大程度上减小拓扑模型的个体差异性的影响。

(4)提出将对一个血管点对进行重建的反投影误差作为它们的势能，通过搜索势能的最大下降方法修正匹配点对，直到达到一个势能稳定状态，可以获得局部最优的匹配结果。

表1是本发明和中国专利文献CN101301207A中的血管匹配准确度对比，其中A，B，...K分别为从四个病人的X射线造影图中提取的11对造影图，A，B，C，D(分别对应病人1、2、3、4)中从左右造影图中提取的血管树拓扑结构相同，E，F，G，H(分别对应病人1、2、3、4)中从左右造影图中提取的血管树拓扑结构不相同，I，J，K(分别对应病人1、2、3)中左右造影图不在心脏运动周期中的同一时刻，x/y表示血管树的总血管段数为y，标记正确血管段数为x；

造影图对	CN101301207A匹配效果	本发明匹配效果
			A	100％	100％
B	46.7％	100％
			C	69.2％	100％
D	76.9％	100％
E	86.7％	100％
			F	46.7％	100％
G	60％	77.8％
			H	66.7％	90.9％
			I	73.3％	100％

造影图对	CN101301207A匹配效果	本发明匹配效果
			J	46.7％	100％
K	69.2％	81.8％

表1

附图说明

图1是X射线造影图血管匹配的流程框图；

图2是X射线造影图及其骨架提取结果；

图2(a)是LCA左视角(-26.8，-27.2)造影图像；

图2(b)是LCA右视角(50.8，30.2)造影图像；

图2(c)是LCA左视角血管分割并提取中轴的结果图；

图2(d)是LCA右视角血管分割并提取中轴的结果图；

图3是由最接近于左冠脉血管树根节点的三个分叉点P1、P2、P3寻找两个血管主枝分叉点的示意图；

图4是模型拓扑结构与造影图的对比图；

图4(a)是模型的拓扑结构；

图4(b)和图4(c)为两种典型对应情况，其中虚线显示的结构都与模型图4(a)的拓扑结构相似；

图5是模型投影图与造影图的相似度比较，其中j点分别为模型投影图和造影图中一个对应分叉点，1和2分别为此分叉点的两个子血管段；

图5(a)为模型投影图；

图5(b)为造影图；

图6是由模型血管树投影图指导造影图血管标记示意图；

图6(a)是指导重建的血管树模型中血管分支的标记示意图；

图6(b)是用模型指导对造影图像中血管分支作标记的结果示意图；

图6(a)中的G1，G2，...，G7分别与图6(b)中的H1，H2，...H7标记相同；

图7是左右造影图中未标记血管的三种对应情况；

图7(a)表示左造影图中未命名血管段UnName1在右造影图中没有对应血管段；

图7(b)表示左造影图中未命名血管段UnName1在右造影图中恰好有一个对应血管段UnName2；

图7(c)表示左造影图中未命名血管段UnName1在右造影图中相同位置处有两个对应血管段UnName2和UnName2’；

图8是计算匹配代价函数时的血管对应情况；

图9是分别对两幅不同视角的血管造影图中的血管进行标记和匹配的最后结果，其中造影图11(a)中的血管段J1，J2，...，J13分别与22(b)中的血管段K1，K2，...，K13匹配对应；

图10是造影图血管点匹配的算法流程图；

具体实施方式

以下结合附图和具体实施例对本发明作进一步说明。

如图1所示，本发明方法利用冠状动脉的多尺度模型对两个不同角度的造影图中的血管段进行迭代匹配，并通过势能下降在对应血管段上获取局部最优的匹配血管点对。具体步骤如下(以左冠脉为例，右冠脉的匹配与此类似)：

(1)X射线造影图中血管拓扑结构提取

通过分割和细化从造影图中提取得到冠脉血管树的骨架，从冠脉血管树的根节点开始，得到血管的树型结构表示，去除血管交叉和伪血管。

图2(a)(b)分别是两个不同角度的造影图，(c)(d)分别是其骨架提取结果。

(2)造影图血管段匹配

根据冠状动脉的解剖学知识，人体的左冠状动脉为一短干，发自左主动脉窦，分为前室间支(也叫左前降支)和左旋支两个主支。首先通过解剖学知识识别和匹配左冠脉的两个主要分支，然后对左右造影图中的这两个主支上的分支进行迭代识别和匹配。

主要包括下述三个步骤：

(2.1)通过造影图中血管的半径信息分别确定左右造影图中左冠状动脉血管树的两个主子树

两个主子树分别为左冠状动脉两个血管主支即前室间支和左旋支所在的血管子树。一般情况下，前室间支和左旋支为左冠状动脉血管树根节点分出的两个血管分支，但是由于左冠状动脉主支上上可能会连接一些小血管，也可能由于噪声影响，经过血管分割得到的冠脉血管树的根节点的第一个分叉点分成的两个血管分支并不是前室间支和左旋支，需要通过血管的半径信息来确定。合理假设最接近于左冠脉血管树根节点的三个分叉点可能为两个血管主枝的分叉点，如图3所示。此血管树中两个血管主枝的分叉点可能为P1或P2或P3，并对其可能性设定不同的优先级，离根节点越近，优先级越高。

令血管树根节点的血管的平均半径为r，根据下列代价函数分别计算点P1，P2和P3作为前室间支和左旋支分叉点的代价。

f_i＝a_i*(|r-lR_i|+|r-rR_i|)

其中a_i为点P_i作为分叉点的优先级系数(i＝1，2或3)，这里设定a₁＝0.6，a₂＝1，a₃＝1，lR_i和rR_i分别为分叉点P_i的左右两个分支。

取代价函数f最小的分叉点为前室间支和左旋支的分叉点。此分叉点的两个分支即为左冠脉血管树的两个主子树，为避免在二维空间中识别错误，暂时先不判断哪个分支是前室降支，哪个是旋支。

(2.2)按X射线造影方式从不同角度对血管树静态模型进行投影，通过投影图与造影图之间的匹配对左右两幅造影图中部分血管进行命名，从而实现左右造影图中部分血管段的匹配。

对每个角度的投影图，通过血管的拓扑结构与位置关系将模型投影图的两个主子树分别与造影图中的两个主子树进行匹配，从而对造影图中的部分血管进行标记。

由于造影图的两个主子树与模型投影图的对应主子树之间可能有多种对应关系，如图4所示列出了两种典型情况，图4(b)和(c)中虚线显示的血管树结构都与(a)中模型表示的血管树相似或相对应，需要在所有对应关系中找到最相似的那种。要搜索所有的对应关系，相当于要在造影图中血管树上找到所有与模型投影图拓扑结构一样的子树。

给造影图中血管树上任意一个分叉点设计一个信号量S_k，约定当S_k＝0时，该分叉点在子树中也是一个分叉点；当S_k＝1时，该分叉点的父血管段与左子血管段在子树中连接成一段血管；当S_k＝2时，该分叉点的父血管段与右子血管段在子树中连接成一段血管。设造影图中血管树上有n个分叉点，对S_k(k＝1，2，...，n)进行3进制编码，则共存在3ⁿ种编码情况，每一种编码情况对应一种可能与模型投影图中血管树拓扑结构相似的血管结构。对每一种编码，判断其对应的血管结构是否与模型投影图中血管树的拓扑结构相似，如果相似即确定了两者是对应的。对每种对应情况，计算左右造影图主子树之间的匹配代价f₁，以及模型投影图主子树与造影图主子树之间的匹配代价f₂。

(2.2.1)计算左右造影图主子树之间的匹配代价f₁：

取f₁为对左右造影图中匹配血管段进行重建后的反投影与原血管段的误差的平均值。

$f_1=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\sqrt{\mathrm{dis}{(P_{1i}-Q_{1i})}^2+\mathrm{dis}{(P_{2i}-Q_{2i})}^2}/N$

其中P_1i和P_2i分别为对应血管段上选取的对应血管点，Q_1i和Q_2i是将P_1i和P_2i进行重建之后按左右造影图参数反投影的点，N是选取的点对的个数。

(2.2.2)计算模型投影图主子树与造影图主子树之间的匹配代价f₂：

f₂＝ω₁*f(angle)+ω₂*f(length)

$f\left(\mathrm{angle}\right)=\underset{i=1}{\overset{M}{\mathrm{\Σ}}}|{\mathrm{angle}}_{1i}-{\mathrm{angle}}_{2i}|$

$f\left(\mathrm{length}\right)=\underset{j=1}{\overset{L}{\mathrm{\Σ}}}\sqrt{|l_{11j}*l_{22j}-l_{12j}*l_{21j}|}$

其中f(angle)表示两者之间血管拓扑结构的匹配代价，angle_1i和angle_2i分别为模型投影图和造影图中对应血管段的第一个点与其中心点连接所成的角度，f(length)表示两段血管长度的匹配代价，l_11j和l_12j，l_22j和l_21j分别为模型投影图和造影图中对应分叉点的左右子血管段的长度。ω₁和ω₂分别为f(angle)和f(length)的匹配代价系数，这里令ω₁＝mean(l_2i)，ω₂＝1，M是模型投影图与造影图对应血管段的数目，L是模型投影图与造影图对应分叉点的数目。

比较每个投影角度下的每种匹配方式的匹配代价f₁和f₂，当f₁具有唯一最小值时，则该f₁所对应的投影角度和匹配方式为最佳投影角度与最优匹配方式，当f₁最小值不唯一时，那么则继续比较在这些情况下的f₂，取f₂最小时对应的投影角度和匹配方式为最佳投影角度与最优匹配方式，然后通过最佳投影角度与最优匹配方式利用模型对造影图中部分血管进行标记。具体算法如下：

初始化最佳投影角度opt_a＝0，opt_b＝0，和最优匹配方式编码opt_code；

初始化最小代价min_f1＝MAX，min_f2＝MAX，MAX为一极大值；

for a＝min_a:step_a:max_a，b＝min_b:step_b:max_b

以角度a和b，以及造影图参数1和D对模型进行投影得到投影图P；

for code＝1:num_code

以编码code对应的方式对模型投影图P和造影图进行匹配，利用模型投影图P中血管名字对造影图中部分血管进行标记，利用匹配结果计算代价函数f₂；利用左右造影图中标记结果对两幅造影图中部分血管进行匹配，利用匹配结果计算代价函数f₁；

if f₁＜min_f1

   min_f1＝f₁

   min_f2＝f₂；

         opt_a＝a；opt_b＝b；opt_code＝code；

    else if f₁＝min_f1

         if f₂＜min_f2

            min_f1＝f₁；

            min_f2＝f₂；

            opt_a＝a；opt_b＝b；opt_code＝code；

         end    

      end

   end

end

图6是血管树模型指导的造影图中血管分支的标记示意图。

(2.3)对于两幅造影图没有标记的血管，根据拓扑结构的约束和最小化反投影误差对它们进行迭代匹配。

设左右两幅造影图中血管拓扑结构分别表示为left_T和right_T，对于所述两幅造影图中主子树上还没有命名标记的所述其余部分血管段，存在图7所示三种对应情况：

对于(a)中情况，left_T中未命名血管段UnName1在right_T中没有对应血管段，则对UnName1不予处理；

对于(b)中情况，left_T中未命名血管段UnName1在right_T中恰好有一个对应血管段UnName2，则令UnName1和UnName2为匹配血管。若UnName1和UnName2不止是一段血管，下面还有分支的话，则需要通过步骤2中的算法搜索得到最优匹配方式；

对于(c)中情况，left_T中未命名血管段UnName1在right_T中相同位置处有两个对应血管段UnName2和UnName2’，则UnName1可能与UnName2匹配，也可能与UnName2’匹配，分别计算这两种可能下的代价函数g，取代价函数最小的为匹配血管。假设血管对应情况如图8所示，代价函数g可通过下式计算：

g＝Err_pro(I1_1)+Err_pro(I1_2)+Err_pro(UnName1，UnName2)

其中，Err_pro(I1_1)和Err_pro(I1_2)表示分别将左右造影图中标记为I1_1和I1_2的血管段进行重建的反投影误差，Err_pro(UnName1，UnName2)表示将UnName1和UnName2进行重建的反投影误差。

图9是分别对两幅不同视角的血管造影图中的血管进行标记和匹配的最后结果，其中造影图9(a)中的血管段J1，J2，...，J13分别与图9(b)中的血管段K1，K2，...，K13匹配对应。

本发明提出的双视角造影图血管段匹配方法，匹配正确率相对于中国专利文献CN101301207A公开的技术方案有了显著的提高。对从四个不同个体的造影图序列中提取的11对造影图分别用本发明方法和上述专利文献中的方法对造影图中的血管进行识别标记，统计标记结果如表1所示。可以看出，无论从左右两幅造影图中分割得到的血管拓扑结构相同(A、B、C和D)还是不同(E、F、G和H)，左右两幅造影图处于心脏运动周期中的同一时刻(A、B、C和D以及E、F、G和H)还是不同时刻(I、J和K)，本文算法都能对造影图中冠状动脉血管进行准确的标记。总的来说，本发明的平均血管匹配正确率能达到96.5％，上述专利文献的平均血管匹配正确率为67.4％。

(3)造影图血管点匹配

本发明提出了一种向最小势能梯度方向迭代匹配血管点的方法，步骤如下：

(1)首先通过外极线约束得到一系列初始的血管匹配点对(P′_i，Q′_i)(i＝1，2，...，m，m为血管匹配点对的个数)；

(2)由于匹配点在血管段上的空间位置需相对应，基于此顺序约束，调整初始匹配点对，得到调整后的血管匹配点对(P_i，Q_i)；

(3)找到调整后的血管匹配点对(P_i，Q_i)中势能梯度最大的点对，在保持所述顺序约束条件下，沿着血管以像素为单位向势能降低最大的方向调整此匹配点对的位置；

其中，点对(P_i，Q_i)的势能即为对这两个点进行三维重建后的反投影误差Err(P_i，Q_i)，而点对(P_i，Q_i)的势能梯度定义为：

Grad(P_i，Q_i)＝max((Err(P_i，Q_i)-Err(P_i，Q_i_pre))，(Err(P_i，Q_i)-Err(P_i，Q_i_next))，

(Err(P_i，Q_i)-Err(P_i_pre，Q_i))，(Err(P_i，Q_i)-Err(P_i_next，Q_i)))

上式中，P_i_pre，P_i_next，Q_i_pre和Q_i_next分别为点P_i和Q_i像素邻域的前点与后点，且它们均在血管段上；

(4)重复上述步骤(3)，直到所有匹配点对都到达一个最低势能状态。

图10为血管点迭代匹配方法的流程图。

Claims (5)

1.一种用于冠脉造影三维重建中的匹配方法，包括血管段的匹配和血管点的匹配，具体步骤如下：

(一)、造影图血管段的匹配，主要包括：

(1)提取基础数据：从两幅不同视角的X射线血管造影图即左右造影图中，分别提取血管中轴图像，从而获得各自的血管树即血管拓扑结构，及血管的半径信息；

(2)预匹配：通过血管半径信息识别并匹配左右造影图中左冠脉或右冠脉的两个主要分支，确定所述两个主要分支各自的主子树，从而确定所述主子树在左右造影图中的对应关系；

(3)对所述主子树进行迭代识别和匹配，具体包括以下两个步骤：

A.模型指导下的部分匹配步骤，即按X射线造影方式从不同角度对血管树静态模型进行投影，通过投影图与造影图之间的匹配对左右两幅造影图中部分血管段进行分别命名标记，从而实现左右造影图中部分血管段的匹配；

B.主子树其余部分的匹配步骤，根据拓扑结构的约束和最小化反投影误差对它们进行迭代匹配，具体匹配过程如下：

设左右两幅造影图中血管拓扑结构分别表示为left_T和right_T，对于所述两幅造影图中主子树上还没有命名标记的所述其余部分血管段，存在三种对应情况：

第一种情况，left_T中未命名血管段UnName1在right_T中没有对应血管段，则对UnName1不予处理；

第二种情况，left_T中未命名血管段UnName1在right_T中恰好有一个对应血管段UnName2，则令UnName1和UnName2为匹配血管，若UnName1和UnName2不止是一段血管，下面还有分支，则通过搜索left_T与right_T中UnName1和UnName2以及它们的分支之间的所有对应关系，利用最小化反投影误差得到最优匹配方式；

第三种情况，left_T中未命名血管段UnName1在right_T中有两个以上对应的血管段，那么未命名血管段UnName1有可能与其中任何一个血管段匹配，因此需分别计算所有可能情况下的代价函数g，取其中代价函数最小的为匹配血管段，代价函数g通过下式计算：

g＝Err_pro(I1_1)+Err_pro(I1_2)+Err_pro(UnName1，UnName2)

其中，UnName2为所述两个以上对应的血管段中任一个血管段，I1_1和I1_2表示left T或right_T中标记为I1_1和I1_2的两个首尾相连的血管段，在left_T中，I1_1和I1_2被UnName1所分，在right T中，I1_1和I1_2被UnName2所分，Err_pro(I1_1)和Err_pro(I1_2)表示分别将左右造影图中标记为I1_1和I1_2的血管段进行重建的反投影误差，Err_pro(UnName1，UnName2)表示将UnName1和UnName2进行重建的反投影误差。

(二)、造影图血管点的匹配，即对匹配好的血管段中的各血管点进行匹配，包括以下步骤：

(1)首先通过外极线约束得到一系列初始的血管匹配点对(P′，Q′)(i＝1，2，...，m，m为血管匹配点对的个数)；

(2)由于匹配点在血管段上的空间位置需相对应，基于此顺序约束，调整初始匹配点对，得到调整后的血管匹配点对(P_i，Q_i)；

(3)找到调整后的血管匹配点对(P_i，Q_i)中势能梯度最大的点对，在保持所述顺序约束条件下，沿着血管以像素为单位向势能降低最大的方向调整此匹配点对的位置；

其中，匹配点对(P_i，Q_i)的势能即为对这两个点进行三维重建后的反投影误差Err(P_i，Q_i)，而匹配点对(P_i，Q_i)的势能梯度定义为：

Grad(P_i，Q_i)＝max((Err(P_i，Q_i)-Err(P_i，Q_i_pre))，(Err(P_i，Q_i)-Err(P_i，Q_i_next))，

(Err(P_i，Q_i)-Err(P_i_pre，Q_i))，(Err(P_i，Q_i)-Err(P_i_next，Q_i)))

上式中，P_i_pre，P_i_next，Q_i_pre和Q_i_next分别为点P_i和Q_i的邻域的前点与后点，且它们均在血管段上；

(4)重复上述步骤(3)，直到所有匹配点对都达到一个最低势能状态。

2.根据权利要求1所述的用于冠脉造影三维重建中的匹配方法，其特征在于，所述的模型指导下的部分匹配步骤具体过程为：

搜索左右造影图的两个主子树与模型投影图的相对应的主子树之间所有对应关系，对每种对应关系，计算左右造影图主子树之间的匹配代价f₁，以及模型投影图主子树与造影图主子树之间的匹配代价f₂，利用最小化匹配代价原则找到最佳投影角度与最优匹配方式，即在造影图的血管树上找到与模型投影图拓扑结构最相似的主子树，进而对造影图中部分血管段进行标记，实现部分匹配。

3.根据权利要求2所述的用于冠脉造影三维重建中的匹配方法，其特征在于，

(I)左右造影图主子树之间的匹配代价f₁的计算为：

取f₁为对左右造影图中匹配血管段进行重建后的反投影与原血管段的误差的平均值：

$f_1=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\sqrt{\mathrm{dis}{(P_{1i}-Q_{1i})}^2+\mathrm{dis}{(P_{2i}-Q_{2i})}^2}/N$

其中P_1i和P_2i分别为对应血管段上选取的对应血管点，Q_1i和Q_2i是将P_1i和P_2i进行重建之后按左右造影图参数反投影的点，N是选取的点对的个数。

(II)模型投影图主子树与造影图主子树之间的匹配代价f₂的计算为：

f₂＝ω₁*f(angle)+ω₂*f(length)

$f\left(\mathrm{angle}\right)=\underset{i=1}{\overset{M}{\mathrm{\Σ}}}|{\mathrm{angle}}_{1i}-{\mathrm{angle}}_{2i}|$

$f\left(\mathrm{length}\right)=\underset{j=1}{\overset{L}{\mathrm{\Σ}}}\sqrt{|l_{11j}*l_{22j}-l_{12j}*l_{21j}|}$

其中f(angle)表示两者之间血管拓扑结构的匹配代价，angle_1i和angle_2i分别为模型投影图和造影图中对应血管段的第一个点与其中心点连接所成的角度，f(length)表示两段血管长度的匹配代价，l_11j和l_12j，l_22j和l_21j分别为模型投影图和造影图中对应分叉点的左右子血管段的长度，ω₁和ω₂分别为f(angle)和f(length)的匹配代价系数，这里令ω₁＝mean(l_2j)，ω₂＝1，M是模型投影图与造影图对应血管段的数目，L是模型投影图与造影图对应分叉点的数目。

4.根据权利要求2或3所述的用于冠脉造影三维重建中的匹配方法，其特征在于，所述的最小化匹配代价原则具体为：比较每个投影角度下的每种匹配方式的匹配代价f₁和f₂，f₁为左右造影图之间的匹配代价，f₂为模型投影图与造影图之间的匹配代价，当f₁具有唯一最小值时，则该f₁所对应的投影角度和匹配方式为最佳投影角度与最优匹配方式，当f₁最小值不唯一时，那么则继续比较在这些情况下的f₂，取f₂最小时对应的投影角度和匹配方式为最佳投影角度与最优匹配方式，然后通过最佳投影角度与最优匹配方式利用模型对造影图中所述的部分血管段进行标记。

5.根据权利要求2或3或4所述的用于冠脉造影三维重建中的匹配方法，其特征在于，所述搜索对应关系的具体过程为：给造影图中血管树上任意一个分叉点设计一个信号量S_k，约定当S_k＝0时，该分叉点在子树中也是一个分叉点；当S_k＝1时，该分叉点的父血管段与左子血管段在子树中连接成一段血管；当S_k＝2时，该分叉点的父血管段与右子血管段在子树中连接成一段血管，设造影图中血管树上有n个分叉点，对S_k(k＝1，2，...，n)进行3进制编码，则共存在3ⁿ种编码情况，每一种编码情况对应一种可能与模型投影图中血管树拓扑结构相似的血管结构，对每一种编码，判断其对应的血管结构是否与模型投影图中血管树的拓扑结构相似，如果相似即确定了两者是对应的。

Images