CN101799935B - 一种单臂x射线造影图像的动态三维重建方法 - Google Patents

Abstract

一种单臂X射线造影图像的动态三维重建方法，属于数字图像处理与医学成像的交叉领域，目的是满足临床医学上心血管疾病的辅助检测以及手术导航的要求。本发明提出了“动态的血管三维重建”概念，对不同时刻的双视角造影图进行呼吸运动补偿和心脏运动补偿。本发明可以获得较高的冠状动脉血管三维重建精度，解决由不同时相的多视角血管造影图进行可靠的自动三维重建难题，有效地辅助心血管疾病的检测与手术导航，满足临床的要求。

Description

一种单臂X射线造影图像的动态三维重建方法

技术领域

本发明属于数字图像处理与医学成像的交叉领域，具体涉及一种单臂X射线造影图像的三维重建方法。

背景技术

血管树三维重建是通过不同视角的X射线二维投影图像中相应的图像信息恢复血管三维空间结构的过程。它与普通可见光图像三维重建有很大不同。X射线成像是投射到人体的X射线经过人体中不同组织的衰减后在荧光屏上形成的图像，每个像素点的值是由在X射线路径上所有组织的衰减的叠加，背景复杂，而且噪声很强，要从X射线造影图重建得到三维心血管树难度很大。

常见的X射线造影成像系统有单臂造影系统和双臂造影系统两种。由于双臂系统保证了两个不同角度的造影图对在时间上的对应性，适合于进行冠脉血管的三维重建，以往的冠状动脉三维重建方法大都是基于双臂造影系统的重建。但由于双臂造影系统占用空间大、成本高，目前临床上大都采用单臂造影系统，通过旋转造影臂得到对应于不同造影角度的造影图或造影图序列。单臂造影系统可以很方便对病人进行不同角度的造影，但是对于重建来说，它存在一个很大的缺陷，即无法得到对应于同一时刻的不同视角的造影图。由于心脏不停地周期运动，而且还是非刚体运动，加上人体的呼吸运动，不同时刻的造影图像中心脏和冠状动脉的形态都发生了变化，还有造影剂到达时刻不同。在此情况下，要获得较高的冠状动脉血管三维重建精度将会十分困难。

为了解决单臂造影重建难题，我们提出了“动态重建”的概念。所谓动态的血管三维重建，就是指被选择用来做血管三维重建的两幅造影图A和B对应于不同时刻。由于冠状动脉血管的运动，造影图A和B对应的冠状动脉的空间三维形态会发生变化，实际上对应两个不同的空间三维形态，此时对造影图A和B的三维重建是一个不适定的问题，目前的静态的血管三维重建方法就不适用了，需要发展新的动态的血管三维重建途径。

目前国内外还没有其他人明确地提出动态的血管三维重建这样的概念，但是也有一些关于血管运动补偿的研究，主要包括两个方面：(1)通过二维图像序列对序列中的某帧图像进行运动补偿。这类方法通常是通过光流或滤波的方法从造影图序列提取血管的二维运动。Curwen利用卡尔曼滤波的方法结合蛇形(snake)模型进行血管运动的估计，Jolly通过光流计算的方法对造影图序列中血管的运动进行估计，Rong将感兴趣血管段用多线段表示，利用纹理信息在相邻造影图进行图形搜索，估计血管二维运动信息并对其进行跟踪补偿。(2)利用从造影图中提取的三维的心脏与呼吸运动模型进行补偿。Shechter等人提出了一种方法从双臂的X射线造影图序列中获得冠脉血管三维的运动，然后通过B样条展开和最小二乘拟合的方法对冠脉三维运动进行分解得到独立的心脏和呼吸运动的模型，并将模型用于同一个病人的其他角度的造影图序列的心脏和呼吸运动补偿。

上述补偿算法对序列图像进行处理然后估计得到血管的二维和三维运动，都是从图像到图像的单纯的图像处理，并没有从冠状动脉血管的三维运动模型出发来指导冠脉的动态的血管三维重建。从图像来谈图像，并没有从本质上解决动态的血管三维重建问题。同时，通过二维图像序列进行运动补偿的方法对成像序列的质量要求较高，这在实际中难以保证，而且也没有考虑到呼吸运动补偿的问题；而利用从造影图中提取三维心脏与呼吸运动模型进行补偿的方法虽然考虑了两种补偿方式，但是需要首先对病人进行双臂造影，得到心脏与呼吸运动的模型，然后利用它来指导单臂造影的重建。由于X射线造影本身是一种介入性检测手段，在实际情况中要尽量避免多次造影，而且这里还要分别用到双臂和单臂两种造影系统，因此在临床上很难得到应用。

本发明是对前中国专利“动态模型指导下的血管造影三维重建方法”(ZL200810047853.5)的深化和扩展。上述专利仅考虑到心脏运动补偿，补偿方法简单且过于依赖三维冠脉运动模型。其做法是直接以模型投影图血管点的运动矢量作为造影图对应点的补偿矢量，因而没有顾及到不同个体心脏运动的特异性。本发明所提出的单臂X射线造影图像的动态重建方法，先将X射线造影图中冠脉血管的运动分为心脏周期运动引发部分与呼吸运动引发部分，然后分别进行补偿，使之更符合实际情况；其过程中，冠脉三维运动模型只提供初始补偿矢量，并在此基础上以心包约束和重建反投影误差构成代价函数，局部搜索找到最优补偿矢量，这样便既利用了模型作为指导，又兼顾了个体心脏运动特异性，提高了重建精度，可以解决由不同时相的多视角血管造影图进行可靠的自动三维重建难题，满足临床医学心血管疾病辅助检测和手术导航的应用要求。

发明内容

本发明提供了一种单臂X射线造影图像的动态三维重建方法，其步骤包括：

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

(2)建立两个不同视角造影图(即左右造影图)血管段间的匹配；

(3)对不同时相的所述左右造影图进行运动补偿，包括呼吸运动补偿和心脏运动补偿；

(3.1)呼吸运动补偿：提取左右两组造影图序列的呼吸运动曲线，确定呼吸运动的呼气末期或吸气末期为补偿参考点，将准备重建的所述左右造影图分别补偿到各自的补偿参考点；

(3.2)按照下列步骤完成左造影图或右造影图的心脏运动补偿：

(3.2.1)建立动态的冠状动脉模型和心脏模型；

(3.2.2)通过点匹配的方式将模型坐标系变换到造影系统坐标系；

(3.2.3)在冠状动脉模型的动态特性视图中选择对应于左造影图或右造影图的造影角度的特性视图，将所述对应的特性视图中的血管段与左造影图或右造影图中的血管段进行匹配，使之相互对应；

(3.2.4)将所述对应的特性视图中对应分叉点所连接的三个血管段上所有点的运动矢量进行平均，得到一个平均运动方向与补偿大小，作为左造影图或右造影图中对应血管分叉点的初始补偿方向与初始补偿大小；

(3.2.5)在确定的搜索区域内通过代价函数寻找最优补偿矢量D_j(j＝1，2，...，n，n为血管分叉点的个数)。其中，代价函数h如下式所示：

$h=\mathrm{Err}\_\mathrm{pro}(P_j\left(t_{c1}\right)+\stackrel{\→}{d})+\mathrm{dis}(C_j,\mathrm{car}\left(t_{c2}\right))$

上式中，

为平移矢量，

为血管分叉点P_i(t_c1)经过心脏缩放和平移后与右造影图或左造影图中匹配的血管分叉点Q_j(t_c2)进行重建的反投影误差，C_j为血管分叉点P_j(t_c1)经过补偿后与血管分叉点Q_j(t_c2)重建的三维点，dis(C_j，car(t_c2))为三维点C_j到t_c2时刻(或时相)的模型心包car(t_c2)的距离。

当平移矢量

使得代价函数h最小时，此时的

即为将血管分叉点P_j(t_c1)补偿到目标时相t_c2(即右造影图或左造影图所处的心脏运动时相)需要的最优补偿矢量D_j。

(3.2.6)根据左造影图或右造影图中血管分叉点的所述最优补偿矢量D_j计算血管段上血管点P′_i(t_c1)的补偿矢量d_i(i＝1，2，...，m，m为血管段上血管点的个数)。

通过下列式子计算d_i：

当血管点P′_i(t_c1)所在血管段上面没有分支时，d_i＝D_k；

当血管点P′_i(t_c1)所在血管段上下均有分支时，d_i＝α*D_k+(1-α)*D_l；

当血管点P′_i(t_c1)所在血管段下面没有分支时，d_i＝D_l。

其中D_k为血管点P′_i(t_c1)所在血管段与下面的分支的交点的补偿矢量，D_l血管点P′_i(t_c1)所在血管段与上面的分支的交点的补偿矢量，l、k为从1到n中取值的整数，α＝l_2start/l_seg，l_2start为血管点P′_i(t_c1)到所在血管段起始点的距离，l_seg为P′_i(t_c1)所在血管段的长度。

令P″_i(t_c2)＝P′_i(t_c1)+d_i，则P″_i(t_c2)即为补偿的结果。

(4)血管点匹配，对补偿后的左右造影图中，拓扑结构匹配的血管段上的各像素点进行匹配，得到匹配像素点。

(5)三维重建，利用步骤(4)中得出的匹配像素点，计算出对应的空间点的三维坐标，即完成三维重建。

进一步地，所述搜索区域通过以下方式确定：设所述血管点分叉点P(t_c1)的初始补偿方向为α，初始补偿大小为l，则补偿的搜索区域为

{c(ρ，θ)|θ∈[α-π/4，α+π/4]，ρ∈[0，3l]}

其中c(ρ，θ)为补偿向量的极坐标表示，ρ为补偿大小，θ为补偿方向，t_c1为补偿初始时相，即左造影图或右造影图所处的心脏运动时相。

本发明方法充分考虑了心脏运动的复杂性对单臂X射线造影图像三维重建的影响，提出了新的动态重建方法。该发明利用傅立叶级数展开的方法从成像序列中分离提取呼吸运动曲线进行呼吸运动补偿，利用冠脉动态模型指导心脏运动补偿。实验证明，使用补偿后的造影图做冠脉三维重建大大提高了重建精度。具体而言，本发明具有如下几个方面的技术效果：

(1)从理论和模型上说明了利用单臂和双臂造影图像重建的差别，论证了利用单臂造影图像进行冠脉血管三维重建的局限性和不完备性，提出了“动态的血管三维重建”的概念，建立了基于单臂造影的重建模型；

(2)动态重建采用呼吸运动补偿。呼吸运动对重建的影响是以往很多方法都忽略的问题，而事实上，对于单臂造影的重建来说，呼吸运动经常会造成很大的重建误差，因此必须进行去除。其他有些方法也考虑了呼吸运动，但是要么通过在人体内植入一些标记点记录人体的呼吸运动，要么首先在双臂造影系统的条件下，分离得到三维的心脏与呼吸运动，然后用于单臂的动态重建，这对单臂造影图像三维重建没有意义。我们提出的动态重建方法是从单臂造影图序列中提取呼吸运动曲线进行呼吸运动补偿，克服了单臂造影系统的局限同时完成了补偿，具有重大意义。

(3)利用模型进行心脏运动补偿。以往的心脏运动估计与补偿都是基于从图像到图像的单纯的图像处理，并没有从冠状动脉血管的三维运动模型出发来指导冠脉的动态的血管三维重建。从图像来谈图像，并没有从本质上解决动态的血管三维重建问题。本文通过心脏与冠脉的三维运动模型对造影图中心脏的运动进行补偿，解决了动态的血管三维重建问题，而且并具有很好的鲁棒性。

附图说明

图1是动态重建方法的流程框图；

图2(a)～图2(d)是原始造影图及中轴提取结果；其中，图2(a)是LCA左视角(-26.8，-27.2)造影图像；图2(b)是LCA右视角(50.8，30.2)造影图像；图2(c)是LCA左视角血管分割并提取中轴的结果图；图2(d)是LCA右视角血管分割并提取中轴的结果图；

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

图4是考虑心脏运动的单臂造影不同时相的两个视角的投影示意图；

图5是病人A在两个不同角度的造影图序列(部分)；其中图5(a)是造影角度为(-26.8°，-27.2°)的序列，图5(b)是造影角度为(50.8°，30.2°)的序列；

图6是呼吸运动曲线，曲线1为横向坐标(X轴坐标)的变化，曲线2为纵向坐标(Y轴坐标)的变化；其中图6(a)是从图5(a)中提取到的呼吸运动曲线，图6(b)是从图5(b)中提取到的呼吸运动曲线；

图7是心脏运动补偿矢量搜索范围示意图；

图8是心脏运动补偿效果；其中，图8(a)为选取的两个不同角度的3幅造影图。I1和I2为从病人A的角度为(-26.8°，-27.2°)的造影图序列中找到心脏舒张末期与收缩末期两个时刻的两幅造影图像，I3为在角度为(50.8°，30.2°)的造影图序列中找到心脏收缩末期的造影图像；图8(b)为将I1进行运动补偿后与I3进行三维重建的反投影对比，其中实线为重建前原图，虚线为反投影图；图8(c)为造影图骨架补偿效果图，其中细实线为原图I1，虚线为补偿后的血管中轴线，粗黑实线为I2；

图9是冠脉造影系统几何模型；

图10是对造影角度为(-26.8°，-27.2°)和(50.8°，30.2°)下的二维呼吸运动进行重建后获得的三维结果；

具体实施方式

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

本发明利用傅里叶级数展开方法提取呼吸运动曲线对二维冠脉造影图进行呼吸运动补偿，利用通过冠脉运动模型构建的运动特性试图指导心脏运动补偿，使两幅单臂X射线造影图像处于相同的心脏与呼吸运动时相，从而进行三维重建。图1是本发明方法的流程框图，具体步骤如下：

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

对X射线造影图进行分割得到冠脉血管的骨架，用八连通链码来进行心血管骨架跟踪，提取血管半径，用二叉树存储血管树的拓扑结构。

图2(a)是LCA左视角(-26.8，-27.2)造影图像；图2(b)是LCA右视角(50.8，30.2)造影图像；图2(c)是LCA左视角血管分割并提取中轴的结果图；图2(d)是LCA右视角血管分割并提取中轴的结果图。

(2)血管匹配

现有技术中存在多种血管匹配方法，如基于知识指导的血管匹配方法，基于拓扑相似性的血管匹配方法，基于模型指导的血管匹配方法以及基于蛇形的血管匹配方法等。在本实施例中所采用的基于模型指导的血管匹配方法。

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

(3)运动补偿

若选取的两幅不同视角的血管造影图A和B对应于不同时刻，需要进行动态的血管三维重建。这里的时刻不仅包括心脏运动时相，还包括呼吸运动时相。在单臂造影系统中，心脏运动与呼吸运动的时相可能都不一样，需要分别进行呼吸与心脏运动补偿从而使两幅造影处于相同的心脏与呼吸运动时相。假设造影图A和B分别对应时刻(t_c1，t_r1)和(t_c2，t_r2)(其中(t_cx，t_rx)表示对应于心脏运动时相t_cx和呼吸运动时相t_rx)，空间点C在(t_c1，t_r1)和(t_c2，t_r2)时刻在造影图A和B上的投影点分别为P(实际上为P(t_c1，t_r1))和Q(实际上为Q(t_c2，t_r2))，如图4所示(图中为方便表示，t₁＝(t_c1，t_r1)，t₂＝(t_c2，t_r2))，因此点P和Q为匹配点对。但是由于它们不对应于同一时刻，它们对应的空间三维点发生了变化，因此不能直接对点P和Q进行重建，需要进行运动补偿。可以将点P(t_c1，t_r1)补偿到t₂时的点P(t_c2，t_r2)，再将点P(t_c2，t_r2)与点Q(t_c2，t_r2)进行重建，得到点P(t_c2，t_r2)与点Q(t_c2，t_r2)对应的空间三维点点C(t_c2，t_r2)，或者将点Q(t_c2，t_r2)补偿到t₁时的点Q(t_c1，t_r1)，再将点Q(t_c1，t_r1)与点P(t_c1，t_r1)进行重建，得到点Q(t_c1，t_r1)与点P(t_c1，t_r1)对应的空间三维点C(t_c1，t_r1)。

下面分两步对血管的运动进行补偿：

(3.1)呼吸运动补偿

呼吸运动补偿的目的就是将两幅不同视角的血管造影图A和B对应的呼吸运动时相t_r1和t_r2变到同一个时相。在三维空间中，呼吸运动使心脏发生三维的平移，在X射线造影图中，呼吸运动对冠状动脉血管的作用就是使其发生二维平移，因此对于呼吸运动的补偿就是对平移的补偿。步骤包括：

(3.1.1)提取呼吸运动曲线

合理假设呼吸运动引起的心脏的运动和与它相邻的器官的运动在造影图平面上是一致的，可以在造影图中找到心脏外的其他组织上的一些特征点作为标记点。在整个序列中跟踪这些标记点，得到这些标记点的运动情况，然后将这些标记点的运动近似为此二维投影面上的呼吸运动。对两个不同视角下的造影图图像序列中的非心脏结构特征点分别进行序列跟踪，获取两个不同视角下的呼吸运动曲线分别记为curve1(t_r)和curve2(t_r)，t_r为呼吸运动的时相。

图5为病人A在两个不同角度的造影图序列(部分)。(a)序列的造影角度为(-26.8°，-27.2°)，(b)序列的造影角度为(50.8°，30.2°)。图6为提取的呼吸运动曲线。图6(a)为从图5(a)中提取到的呼吸运动曲线，图6(b)为从图5(b)中提取到的呼吸运动曲线。曲线1为横向坐标(X轴坐标)的变化，曲线2为纵向坐标(Y轴坐标)的变化。

(3.1.2)确定补偿参考点

找到呼吸运动的两个顶点(呼气末期和吸气末期)，选取其中一个点作为补偿参考点(根据实际情况决定选谁作为参考点，有些序列里面可能只有吸气末期，或呼气末期，可能的话最好是选择呼气的点)，在将左右两幅造影图都补偿到所述参考点，这样这两幅造影图就可以看作在呼吸运动的同一时刻，从而完成呼吸运动补偿。

步骤(3.1.1)中提取得到的呼吸运动向量随时间的变化曲线分别为curve1(t_r)和curve2(t_r)，在左右造影图序列中找到了同一个呼吸运动参考时相t_r0(即补偿参考点，t_r0对应于呼气末期或吸气末期)。图6中直线标记的为选取的呼吸运动参考时相，图6(a)中选取的参考帧是第57帧，(b)中选取的参考帧是第28帧。通过下面的公式进行呼吸运动补偿：

P′(t_r0)＝P(t_r1)-(curve1(t_r1)-curve1(t_r0))

Q′(t_r0)＝Q(t_r2)-(curve2(t_r2)-curve2(t_r0))

式中，P′(t_r0)、Q′(t_r0)分别为P(t_r1)、Q(t_r2)经呼吸运动补偿后的点。

(3.2)心脏运动补偿

心脏运动补偿的目的就是将两幅造影图对应的心脏运动时相t_c1和t_c2变换到同一个时相，这里选择将左造影图从t_c1时刻补偿变换到t_c2时刻(右造影图的补偿与此类似)。

(3.2.1)建立动态的冠状动脉模型和动态的心脏模型。通过从心脏切片中提取冠状动脉和心脏信息，重建出三维的静态动脉和心脏模型，在此基础上再加入冠状动脉和心脏的运动信息，从而建立动态的冠状动脉模型和心脏模型。

(3.2.2)坐标变换

坐标变换的目的是将冠状动脉模型的数据从模型坐标系变到造影系统坐标系。由于模型坐标系与造影系统坐标系没有任何关系，只能通过点匹配的方式获得坐标变换的旋转平移矩阵。

从造影图序列中选择舒张末期的两幅不同角度的造影图，进行呼吸运动补偿后重建得到其冠状动脉三维结构，在上述冠状动脉三维结构和冠状动脉模型上找到一些特征点(分叉点)作为对应血管点。假设冠状动脉模型中的血管点Mi与上述冠状动脉三维结构中的血管点Pi对应，则有

X_Mi*R+t＝X_Pi

其中X_Mi＝[x_Mi，y_Mi，z_Mi]，为点Mi的坐标，X_Pi＝[x_Pi，y_Pi，z_Pi]为点Pi的坐标，R为3*3的旋转矩阵，t为1*3的平移矩阵。

上式转化为矩阵形式为

$\left[X_{\mathrm{Mi}}\right|1\left]\right[\frac{R}{t}]=X_{\mathrm{Pi}}$

令A＝[X_Mi|1]，B＝X_Pi，若选取了s对对应点，则A为s*4的矩阵，B为s*3的矩阵，通过最小二乘的方法解方程，得到

$\left[\frac{R}{t}\right]={\left(A^{T}A\right)}^{-1}{A}^{T}B$

得到R，t后，通过旋转平移变化将血管模型从模型坐标空间变换到造影系统的坐标空间。

(3.2.3)确定左造影图血管分叉点的最优补偿矢量D_j

1)在冠状动脉模型的动态特性视图中选择对应于左造影图的造影角度的特性视图，将冠脉特性视图中的血管与左造影图中的血管进行匹配，使之相互对应。

2)选取左造影图中冠状动脉血管所有的分叉点；将特性视图中对应分叉点所连接三个血管段上所有点的运动进行平均得到一个平均运动方向与补偿大小，作为左造影图中血管分叉点的初始补偿方向与初始补偿大小；若造影图中的血管在特性视图中不存在对应血管，则将此血管的初始补偿方向与大小设为相邻血管段的初始运动方向与初始补偿大小。

3)确定搜索区域：设血管分点P_j(t_c1)(j＝1，2，...，n，n为分叉点的个数)的初始补偿方向为α，大小为l，则补偿的搜索区域为

{c(ρ，θ)|θ∈[α-π/4，α+π/4]，ρ∈[0，3l]}

其中c(ρ，θ)为补偿向量的极坐标表示，ρ为补偿大小，θ为补偿方向，搜索区域如图7所示。

4)选取代价函数h如下式所示，在区域C内搜索得到使代价函数h最小的矢量为将血管分叉点P_j(t_c1)补偿到t_c2需要的补偿量。

$h=\mathrm{Err}\_\mathrm{pro}(P_j\left(t_{c1}\right)+\stackrel{\→}{d})+\mathrm{dis}(C_j,\mathrm{car}\left(t_{c2}\right))$

其中

为平移矢量，

为血管点分叉点P_j(t_c1)经过心脏缩放和平移后与右造影图中对应的血管分叉点Q_j(t_c2)进行重建的反投影误差，C_j为血管分叉点P_j(t_c1)经过补偿后与血管分叉点Q_j(t_c2)重建的三维点，dis(C_j，car(t_c2))为点C_j到t_c2时刻的心脏模型的心包car(t_c2)的距离。

(3.2.4)确定左造影图血管段的补偿方向和大小

上一步得到造影图所有血管分叉点的最优补偿矢量D_j，通过下列式子计算左造影图中冠脉血管上每一个血管点P′_i(t_c1)的补偿矢量d_i(i＝1，2，...，m，m为血管段上血管点的个数)：

当血管点P′_i(t_c1)所在血管段上面没有分支时，d_i＝D_k；

当血管点P′_i(t_c1)所在血管段上下均有分支时，d_i＝α*D_k+(1-α)*D_l；

当血管点P′_i(t_c1)所在血管段下面没有分支时，d_i＝D_l。

其中D_k为血管点P′_i(t_c1)所在血管段与下面的分支的交点的补偿矢量，D_l为血管点P′_i(t_c1)所在血管段与上面的分支的交点的补偿矢量，l、k为从1到n中取值的整数，α＝l_2start/l_seg，l_2start为血管点P′_i(t_c1)到所在血管段起始点的距离，l_seg为血管点P′_i(t_c1)所在血管段的长度。

令P″_i(t_c2)＝P′_i(t_c1)+d_i，则P″_i(t_c2)即为补偿的结果。图8是心脏运动补偿效果示意图。

(4)血管点匹配，对补偿后的左右造影图中，拓扑结构匹配的血管段上的各像素点进行匹配，得到匹配像素点。

血管点的匹配方法是本领域的现有技术，属于本领域中血管点匹配的一种常规方法，在本实施例中所采用的是在左右造影图中对应的血管段上进行等间隔取样的方法实现血管点匹配的。

(5)三维重建

血管的三维重建就是找到两个不同角度的冠脉造影图像匹配的血管点对利用冠脉造影系统几何模型，通过空间几何与坐标关系，计算出空间点的三维坐标，即利用步骤(4)中得出的匹配像素点，计算出对应的空间点的三维坐标，完成三维重建。图9所示为冠脉造影系统几何模型。

Claims (2)

1.一种单臂X射线造影图像的动态三维重建方法，其步骤包括：

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

(2)建立两个不同视角造影图的血管段间的匹配，其中，所述两个不同视角造影图分别称为左造影图和右造影图；

(3)对不同时相的所述左右造影图进行运动补偿，包括呼吸运动补偿和心脏运动补偿；

(4)血管点匹配，对补偿后的左右造影图中，拓扑结构匹配的血管段上的各像素点进行点匹配，得到匹配像素点；

(5)三维重建，利用步骤(4)中得出的匹配像素点，计算出对应的空间点的三维坐标，即完成三维重建；

步骤(3)中，所述的呼吸运动补偿为：提取左右两组造影图序列的呼吸运动曲线，确定呼吸运动的呼气末期或吸气末期作为补偿参考点，将准备重建的所述左右造影图分别补偿到各自的补偿参考点；

步骤(3)中，所述的心脏运动补偿具体步骤为：

(I)建立动态的冠状动脉模型和心脏模型；

(II)将模型坐标系变换到造影系统坐标系；

(III)在冠状动脉模型的动态特性视图中选择对应于左造影图或右造影图的造影角度的特性视图，将所述对应的特性视图中的血管段与左造影图或右造影图中的血管段进行匹配，使之相互对应；

(IV)将所述对应的特性视图中对应分叉点所连接的三个血管段上所有点的运动矢量进行平均，得到一个平均运动方向与补偿大小，作为左造影图或右造影图中对应血管分叉点的初始补偿方向与初始补偿大小；

(V)在确定的搜索区域内通过代价函数寻找最优补偿矢量D_j(j＝1，2，...，n，n为血管分叉点的个数)，其中，代价函数h如下式所示：

$h=\mathrm{Err}\_\mathrm{pro}(P_j\left(t_{c1}\right)+\stackrel{\→}{d})+\mathrm{dis}(C_j,\mathrm{car}\left(t_{c2}\right))$

上式中，为平移矢量，为血管分叉点P_j(t_c1)经过心脏缩放和平移后与右造影图或左造影图中匹配的血管分叉点Q_j(t_c2)进行重建的反投影误差，C_j为血管分叉点P_j(t_c1)经过补偿后与血管分叉点Q_j(t_c2)重建的三维点，dis(C_j，car(t_c2))为三维点C_j到t_c2时刻的心脏模型的心包car(t_c2)的距离，

当平移矢量使得代价函数h最小时，此时的

即为将血管分叉点P_j(t_c1)补偿到目标时相t_c2需要的最优补偿矢量D_j；

(VI)根据左造影图或右造影图中血管分叉点的所述最优补偿矢量D_j计算血管段上血管点P′_i(t_c1)的补偿矢量d_i(i＝1，2，...，m，m为血管段上血管点的个数)；

通过下列式子计算d_i：

当血管点P′_i(t_c1)所在血管段上面没有分支时，d_i＝D_k；

当血管点P′_i(t_c1)所在血管段上下均有分支时，d_i＝α*D_k+(1-α)*D_l；

当血管点P′_i(t_c1)所在血管段下面没有分支时，d_i＝D_l；

其中D_k为血管点P′_i(t_c1)所在血管段与下面的分支的交点的补偿矢量，D_l为血管点P′_i(t_c1)所在血管段与上面的分支的交点的补偿矢量，l、k为从1到n中取值的整数，α＝l_2start/l_seg，l_2start为血管点P′_i(t_c1)到所在血管段起始点的距离，l_seg为P′_i(t_c1)所在血管段的长度，

令P″_i(t_c2)＝P′_i(t_c1)+d_i，则P″_i(t_c2)即为补偿的结果。

2.根据权利要求1所述的动态三维重建方法，其特征在于，所述的搜索区域通过以下方式确定：

设所述血管分叉点P_j(t_c1)的初始补偿方向为α，初始补偿大小为l，则补偿的搜索区域为

{c(ρ，θ)|θ∈[α-π/4，α+π/4]，ρ∈[0，3l]}

其中c(ρ，θ)为补偿向量的极坐标表示，ρ为补偿大小，θ为补偿方向，t_c1为补偿初始时相，即左造影图或右造影图所处的心脏运动时相。

Images