动态模型指导下的血管造影三维重建方法

摘要

动态模型指导的血管造影三维重建方法，属于数字图像处理与医学成像的交叉领域，目的是满足我国临床医学上心血管疾病的辅助检测以及手术导航的特殊要求。本发明包括造影图预处理步骤，血管分割步骤，血管骨架与半径提取步骤，模型指导血管基元识别，血管匹配以及血管三维重建步骤。本发明还提供了一种心血管动态模型建立方法，包括心血管切片数据提取步骤，心脏的静态与动态建模步骤，以及心血管系统的静态与动态建模步骤。本发明可以得到很好的血管造影三维重建结果，有效地辅助心血管疾病的检测与手术导航，满足临床的要求。

说明书

技术领域

本发明属于数字图像处理与医学成像的交叉领域，具体涉及一种动态模型指导下的血管造影三维重建方法。该方法可以解决由多视角血管造影图进行可靠的自动三维重建难题，满足临床医学心血管疾病辅助检测和手术导航的应用要求。

背景技术

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

目前医院大都用的是X射线单臂造影系统对病人做X射线造影，通过旋转造影臂得到一个对应于不同造影角度的造影图序列。单臂造影可以很方便对病人进行不同角度的造影，但是对重建来说有一个缺点就是我们无法得到同一时刻的不同视角的造影图，这给重建带来很大困难。

要重建血管树的真实三维空间结构，需要得到血管树至少两个不同角度的投影信息。传统方法首先提取出血管的骨架，然后通过不同视角空间约束关系，对不同视角投影图像的血管像素点进行正确匹配并重建，当血管几何形变不大、几何关系明显的时候，才能大体上恢复出血管的三维空间结构。

通常情况下基于两幅单平面造影图像的血管树的三维重建技术由以下3个关键步骤组成：(1)血管骨架的提取；(2)特征点的识别和匹配；(3)血管空间点的估计和重建及血管基元的拟合。

血管树从整体上可以看作在空间中弯曲延展的管状系统，其骨架是具有树状结构的连续空间曲线，反映了血管树的整体形态特征。通过血管轮廓检测，实现血管骨架和半径准确的自动提取是一个困难的问题。目前，这个方面比较可行的方法有Hoffman等人提出的双四边形区域搜索(double square box region of search)方法，用于自动跟踪血管，同时得到血管的局部大小、中轴线位置等参数。Coatrieux和Collorec等人将矢量跟踪算法用于血管轴线及边界的提取，具有人工干预少，能实现血管由粗到细的分层递进检测和节省时间等优点。但这些方法对原始图像质量要求高，由于目前x光造影成像质量有限，无论是清晰度，还是分辨率，利用现有的图像预处理方法对其增强，效果不佳。另外，要想对血管进行准确的检测与匹配还需要大量的解剖学知识，如果引入专家系统将增加问题的复杂性和成本。临床中大多采用人工检测、人机交互的方法来进行血管骨架的提取。

对于心血管双视角造影图像的结构识别和匹配问题，存在以下几个难点：第一，两个造影成像面中的血管拓扑结构往往相差较大；第二，投影图像中血管投影交叉的存在，严重影响着血管拓扑结构的识别和描述。

传统的方法是用外极线约束+线性规划的方法来对两幅造影图中的血管点进行匹配，但是我们知道，外极线利用的是投影点之间的空间几何关系，它假定的是两个投影点对应的是同一个三维空间点。在单臂造影中我们无法得到对应于同一时刻的造影图，而由于心脏的运动血管点也会发生运动，这样用外极线约束找到的匹配点就会有很大误差。到目前为止，对于血管的三维重建中不同视角血管像素点匹配最常用的方法还停留在人工选取匹配特征点对的方法。这不仅需要操作人员具有丰富经验，而且还会花费较长时间。要有效的识别出双视角造影图像中的匹配的特征点对，并建立不同视角血管匹配的快速而有效的自动算法，需要综合考虑血管方向矢量的连续性、血管的直径信息以及血管拓扑结构约束等等因素。

鉴于传统血管结构识别与匹配方法中存在的问题，国外也曾有人提出了用血管树模型指导造影图匹配的方法，其中大部分用的是根据Dodge提出的由37个个体特异的冠状动脉模型建立的一个通用冠状动脉模型Coronix。由于用来指导造影图匹配的血管树模型是静态的，国内外的很多研究者大都通过在序列图像中选择处于心脏运动周期中同一时刻的造影图作为重建的参考图，比如说舒张末期。但是这样存在两个问题：(1)如何对两幅造影图是否处于一个心脏周期中同一个时刻进行判决。通常判断血管向外舒张到最大(这时血管包着的心脏体积最大)的时候为舒张末期，血管收缩到最小(这时血管包着的心脏体积最小)的时候为收缩末期，但是用人眼判断仍然存在不确定性；(2)这样会浪费掉在除了这个选定时刻外的其他时刻的造影图片，可能会丢掉很多有用的信息。而且这样也没有考虑到造影图的效果，比如造影剂到达时刻、血管遮挡、噪声等，很难保证选定时刻的这些不同方向的血管造影图适合于做三维重建。

在找到不同视图中血管像素点的正确对应关系之后，才可以通过空间几何关系法或者是齐次坐标变换法恢复出血管对应像素点的空间三维坐标。从有限的X-射线投影视图恢复出血管像素点的空间坐标，这种方法只能是一种相似性重建，因而往往通过所重建血管结构在所有视图中的投影与其在对应像面中血管结构之间的欧式距离来构建能量函数，通过最优化方法，找到一个最优的值，即为重建结果。

在进行血管三维重建时，一般不需要重建出血管上每一个点的三维坐标，可以先对造影图中的每个血管段上的进行采样得到一些取样点，然后对这些取样点进行三维重建得到一些三维空间点，再根据插值或曲线拟合的方法从这些三维空间点构建得到整段血管。曲线拟合的方法多种多样，可以根据血管的特殊形态构建适当的曲线拟合模型，较准确的描述血管骨架的空间形态。传统方法中，曲线在空间拟合后可将其反投影到像面，通过比较实际视图中曲线的拟合程度，再对空间血管树模型进行修正，但这个过程需要人工干预，而且准确度并不高。现今成熟的算法有B样条、Snake模型等，可以较准确的表达不同血管基元之间的血管骨架信息。对于样条拟合，实际上是在假定曲线的局部平滑性的基础上在取样点之间进行插值，拟合点越多，精度越高；而对于Snake模型，需要构造合适的内力和外力函数，外力的作用是使拟合曲线尽量逼近投影，内力的作用则是使曲线尽量平滑。相对来说，Snake模型能够获得较高的精度，但内力和外力函数构造相对比较困难；样条拟合的方法比较简单，但提高精度需要较大的计算量。

发明内容

本发明的提供一种动态模型指导的血管造影三维重建方法，目的是通过一个心脏与心血管的四维动态模型，指导双视角造影图中的血管基元投影匹配，进行自动化的、鲁棒的三维重建。

本发明提供的一种基于动态模型的血管造影三维重建方法，其步骤包括：

(1)选取两幅不同视角的X-射线冠状动脉造影图，对它们进行预处理，首先对各造影图图像序列采用频域增强的方法增强血管，然后采用形态学Bottom-Hat变换进一步强调血管，消除干扰噪声，保留并突出所需的图像信息；

(2)对预处理后的造影图像进行分割提取血管区域，细化血管区域得到血管骨架，选用八连通链码来进行心血管骨架跟踪，提取血管半径，按照下述过程采用二叉树结构存储血管树的拓扑结构；

(2.1)在血管造影图的血管中轴提取图像中，从血管树根结点开始对血管中轴进行二叉树搜索，并对搜索过的血管像素做标记，当搜索到标记过的血管像素时表明搜索出现循环，判断有血管发生交叉，并停止向前搜索，转到步骤(2.2)；若搜索过程没有循环出现，转到步骤(2.4)；

(2.2)当二叉树搜索出现循环时，根据血管树的拓扑结构信息和几何信息，在血管造影图的血管中轴提取图像中识别出交叉的血管；

(2.3)回到步骤(2.1)；

(2.4)搜索结束；

(3)按下述过程分别对两幅不同视角的血管造影图的血管中轴提取图像中的血管进行标记：

(3.1)按X射线造影方式从不同方向对选取的血管树模型进行投影得到不同视角的模型投影图；

(3.2)按下述过程求每个投影图与血管造影图的血管中轴提取图像T的相似度，找到其中相似度最高的投影图P；

(3.2.1)将血管树模型投影图和血管造影图的血管中轴提取图像T中的血管拓扑结构分别保存在二叉树T1和T2中；

(3.2.2)从二叉树T2中提取与T1拓扑结构匹配的子树T3；对T1和T3进行匹配，利用下式计算下列T1和T3的代价函数：

f＝a₁·f₁(angle)+a₂·f₂(curve)+a₃·f₃(length)

其中f₁(angle)表示两者之间血管拓扑结构的相似性，f₂(curve)表示两者之间血管结构曲线形状的相似性，f₃(length)表示两段血管长度的相似性，表示为血管长度length的函数；a₁，a₂，a₃为这三种相似性的加权系数，a₁+a₂+a₃＝1，取其中代价函数最小的血管对为匹配血管；

(3.2.3)取其中代价函数最小的子树为T1的匹配子树M，将此最小代价函数的倒数作为模型投影图与血管造影图之间的相似度；

(3.3)通过上述求相似度的步骤得到与血管造影图的血管中轴提取图像T的相似度最高的投影图P在T中的匹配子树M；用P中血管基元投影的名字对子树M中的血管基元投影进行标记，当选取的两个造影图不在心动周期中的同一时刻时，分别用这两个时刻的血管树模型对相应造影图的血管基元投影进行标记；

(4)按下述过程进行血管匹配：

(4.1)找到两幅造影图中相同标记的血管基元投影作为匹配血管基元投影；

(4.2)对照两幅造影图，根据拓扑结构的限制在两幅造影图中没有标记的血管基元投影里找到拓扑结构相同的血管基元投影作为匹配的一对血管基元投影；其过程为：

(4.2.1)设这两幅造影图中血管拓扑结构分别存在left_T和right_T中；

(4.2.2)找到left_T中没有标记的血管UnName1，在right_T中相同位置处搜索是否也有未标记的血管结构UnName2，如果有就说明可能存在对应的血管结构，没有则判断UnName1为伪血管；

(4.2.3)若UnName2存在，判断UnName1和UnName2之间是否一一对应；若是一一对应，则直接将UnName1和UnName2命以相同标记；若非一一对应，则需要两两之间进行匹配，找到代价函数最小的血管对为匹配血管，并作相同标记；

(4.2.4)重复步骤(4.2.1)-(4.2.3)，直到血管树里的血管全部标记或判断为伪血管为止；

(5)通过动态模型指导将两幅造影图中提取的血管中轴点坐标变换到同一时刻；对两幅造影图中拓扑结构匹配的血管基元投影对上的各像素点进一步进行匹配，得到匹配像素点，进行三维重建。

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

(1)使用动态模型指导，增强不同投影角度下血管结构匹配的鲁棒性，同时补偿了心脏运动的影响。在以往的三维重建算法中，大都利用外极线的方法来进行血管的匹配，但是这种方法的误匹配率很高，特别当不同角度获取的造影图不是在同一时刻得到时。为了克服外极线匹配方法的不足，有人提出了用血管树模型指导造影图匹配的方法，虽然能部分地提高血管匹配的精度，但是它仍然存在一个问题：即由于用来指导造影图匹配的血管树模型是静态的，用它来对不同时刻的血管造影图进行匹配得到的结果往往会有很大的误差，特别是血管运动造成血管出现重叠和交叉的时候。这里我们使用动态的血管树模型对血管匹配进行指导，通过不同时刻的模型投影图与选取用来重建的血管造影图之间的拓扑结构匹配分别对造影图中的血管进行标记，克服了血管运动带来的误差。

另一方面，由于进行重建的两幅造影图不在同一时刻，也不能用一般的三维重建公式对它们进行重建。因此我们利用了血管动态模型对造影图中血管进行运动补偿，尽量使两幅造影图像中的血管位置对应于同一个时刻，然后再利用血管三维重建公式对其进行重建，这样可以减少血管运动给重建带来的误差。

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

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

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

附图说明

图1是本发明的流程框图；

图2是心血管建模部分的流程框图；

图3是X射线造影重建的流程框图；

图4是建立心脏及血管树模型的坐标系；

图5(a)是心脏切片图像；

图5(b)是在心脏切片图像中提取心脏部分的图像，其中区域A2是左心室，区域A3是右心室，区域A1是左心房，区域A4是右心房；

图6是心脏轮廓取样点图；

图7是心脏三维模型图，其中最外面的轮廓是心包，区域B1是左心室，区域B2是右心室，区域B4是左心房，区域B3是右心房；

图8(a)～8(c)是通过改变控制点的位置模拟具有个体特异性的心脏图；

图8(a)是正常心脏图；

图8(b)是心脏局部扩张图；

图8(c)是心脏局部内缩图；

图9是心动周期各时期图

图10是左心室容积线性变化图；

图11(a)～11(g)心脏在心动周期各时期的模型图；

图11(a)是心脏处于等容收缩期模型图；

图11(b)是心脏处于快速射血期模型图；

图11(c)是心脏处于缓慢射血期模型图；

图11(d)是心脏处于等容充盈期模型图；

图11(e)是心脏处于快速充盈期模型图；

图11(f)是心脏处于缓慢充盈期模型图；

图11(g)是心脏处于心房收缩期模型图；

图12(a)是心脏切片序列中第48个切片图像；

图12(b)是对(a)提取的血管截面图像；

图12(c)是心脏切片序列中第50个切片图像；

图12(d)是对(c)提取的血管截面图像；

图13是经过血管基元标记的模型示意图；

图14(a)～14(g)是冠状动脉血管在心动周期各时期的模型图示例；

图14(a)是左冠状动脉血管处于心脏等容收缩期模型图示例；

图14(b)是左冠状动脉血管处于心脏快速射血期模型图示例；

图14(c)是左冠状动脉血管处于心脏缓慢射血期模型图示例；

图14(d)是左冠状动脉血管处于心脏等容充盈期模型图示例；

图14(e)是左冠状动脉血管处于心脏快速充盈期模型图示例；

图14(f)是左冠状动脉血管处于心脏缓慢充盈期模型图示例；

图14(g)是左冠状动脉血管处于心脏心房收缩期模型图示例；

图15是左冠状动脉血管和背景噪声点灰度在选取的一个造影图序列中的变化趋势，其中曲线1表示坐标为(252，134)的血管像素的灰度值，曲线2表示坐标为(362，95)的背景点的灰度值。从图中可以看出血管点的灰度变化比背景噪声点要大，根据这个特点可将血管和背景区分开；

图16(a)是左冠状动脉血管造影图像；

图16(b)是左冠状动脉血管造影增强图像；

图16(c)是左冠状动脉血管造影Bottom-Hat变化图像；

图17是血管搜索出现循环示意图；

图18是血管投影交叉示意图；

图19是两幅造影图中血管拓扑结构比较示意图，其中若血管段E1，E2分别和F1，F2对应，则判断血管段UE和UF的拓扑结构相同，否则不同；

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

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

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

图21是两个未标记血管在血管树中同一个拓扑位置示意图，其中血管段I1，I2，I3为相同标记，因此判断未标记血管段U1，U2的拓扑结构相同；

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

图23(a)～图23(i)是通过模型指导冠状动脉血管三维重建的过程以及与其它血管重建算法的比较；

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

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

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

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

图23(e)是与LCA左视角拓扑结构匹配的模型投影图；

图23(f)是与LCA右视角拓扑结构匹配的模型投影图；

图23(g)是用我们的方法对(a)，(b)中造影图进行血管三维重建的结果；

图23(h)是使用常规极线约束的方法对(a)，(b)中造影图进行血管三维重建结果，重建失败；

图23(i)是通过手动选取匹配点的方法对(a)，(b)中造影图进行血管造影三维重建结果；

图24(a)是23(g)中血管三维重建结果按图23(a)中造影角度投影的图像；

图24(b)是23(g)中血管三维重建结果按图23(b)中造影角度投影的图像；

图24(c)和(d)与图23(c)和(d)相同；

图25(a)～(c)是血管多尺度模型示意图；

图25(a)是左冠状动脉一级血管树模型示意图；

图25(b)是左冠状动脉二级血管树模型示意图；

图25(c)是左冠状动脉三级血管树模型示意图；

图26(a)是对附图23中的造影图使用本发明方法重建得到血管三维结构；

图26(b)是对图26(a)中血管中轴加上半径的结果；

图26(c)是对根据虚拟人切片数据建立的血管树模型；

图26(d)是用(a)中重建结果对(c)中模型进行补充细节之后的结果，其中用圆圈标识出的血管段为补充到模型中细尺度的血管基元；

图27是造影系统在两个不同角度的成像示意图。

具体实施方式

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

如图1所示，本发明方法是利用心脏与心血管动态模型从双视角X-射线造影图像中进行血管三维重建，所利用的心脏与心血管动态模型可以采用原有的已知心脏与心血管树模型，也可以采用各人自建的心脏与心血管动态模型，下文提供了一种建立心脏与心血管动态模型的方法，以作参考。

(1)建立心血管树模型

如图2所示为建立心血管树模型流程图。

在建立心脏与血管的三维静态和动态模型时，按如下方式建立坐标系：x、y轴及原点与心尖的切片图的x、y轴及原点重合，z轴垂直于切片平面，沿从心尖到心底的方向。(如图4所示)

(1.1)建立心脏动态模型

建立心脏动态模型的目的是为了帮助建立冠状动脉四维动态模型，其步骤包括：

(1.1.1)提取心脏切片轮廓；

根据解剖学的知识和人体横断解剖的图片，在每幅原始图片中分出心包，左心室，右心室，左心房，和右心房五个部分。

原始数据图片来自于美国可视化人计划(VHP：Visible Human Project)数据集。我们选取的图片是Visible Man数据集中的解剖图像的心脏部分。图5(a)是我们得到的原始心脏切片图。

心脏部分的结构比较复杂，我们通过photoshop来进行分割。图5(b)是我们对图5(a)中的心脏切片图的分割结果。最外面的轮廓是心包，区域A2是左心室，区域A3是右心室，区域A1是左心房，区域A4是右心房。

(1.1.2)对上述心脏切片轮廓进行采样，得到采样点。

由于数据太多会增大计算量，而大部分情况下不需要太多的数据来进行曲线拟合，因此首先对原始数据进行采样，得到一些有代表性的样点。考虑到心脏腔室和表面平滑性的特点，在心脏切片轮廓方向(横向)与垂直于切片方向(纵向)两个方向上都进行平均采样。

横向上的均匀采样是首先找出轮廓的中心(轮廓的形心)，然后从中心发射n条射线，其中相邻射线间的夹角均为360/n度，选取这些射线与轮廓的交点为采样点，因此采样点的个数也为n。图6是取n＝12时得到的12个均匀采样点。

由于纵向(z轴)上有138幅切片图，可以通过每隔几个切片取一幅切片图，得到一些取样图片来进行建模。

需要注意的一点是，取样点越少，得到的模型越平滑；取样点越多，保留的细节越多，但计算量也越大。

(1.1.3)利用上述取样点，根据下述公式进行B样条曲面拟合，建立心脏三维静态模型：

B样条曲面拟合公式为$P(u,v)=\underset{i=1}{\overset{N_u}{\Sigma }}\underset{j=1}{\overset{N_v}{\Sigma }}{C}_{\mathrm{ij}}{B}_i\left(u\right)B_j\left(v\right),$P是模型上的点，u，v是对应此P点的横向(U向)和纵向(V向)参数，N_u，N_v分别是B样条模型中U，V方向上的控制点个数。B_i(u)，B_j(v)则分别是在U，V方向上的B样条基函数，C_ij是B样条曲面的控制点，构成N_u×N_v的控制点网格，决定了B样条曲面的形状。

可以将B样条曲面看作一个由两个方向的B样条曲线进行张量积后的曲面。因此可以将曲面拟合过程分为下列两步：

1)u向的B样条曲线拟合；

B样条曲线拟合公式为$P\left(u\right)=\underset{i=1}{\overset{N_u}{\Sigma }}{C}_i{B}_i\left(u\right),$这里的拟合数据为每一个心脏切片轮廓上的取样点。对这些点进行B样条曲线拟合得到切片轮廓曲线的控制点C′_ij(j＝1，2，...，N_u)。

2)v向的B样条曲线拟合。

这一步的拟合数据为上一步得到的切片轮廓曲线的控制点C′_ij(j＝1，2，...，N_v)。对它们进行B样条曲线拟合得到整个B样条曲面的控制点C_ij(i＝1，2，...，N_u，j＝1，2，...，N_v)。

当得到B样条曲面的控制点后，代入B样条曲面方程，就得到心包和各个腔室的参数模型。指定一个参数值(u₀，v₀)，就可以通过曲面拟合公式计算得到其在空间中对应的点P(u₀，v₀)的坐标(x₀，y₀，z₀)。可以通过改变控制点位置改变心脏形状，模拟具有个体特性的心脏。

图7是将心脏腔室与表面合在一起构成心脏的三维静态模型。最外面的轮廓是心包，区域B1是左心室，区域B2是右心室，区域B4是左心房，区域B3是右心房。

图8是通过改变控制点的位置模拟的具有个体特异性的心脏。

(1.1.4)建立心脏四维动态模型：通过建立心脏各腔室和心包的运动模型来建立动态的心脏模型。

心脏的运动是十分复杂的，根据医学观察和生物医学工程的实验结果证明，室壁运动由心肌的内向收缩运动、心脏的水平移动及心脏沿轴的转动三种主要运动方式组成此外还有扭转、局部伸展等局部运动。在不同的运动形式中，作为心脏运输血液的动力之源，膨胀/收缩运动占到了将近90％，因此我们主要考虑心脏的这种运动方式。在建模前我们先做如下假设：

①左右心房和心室的运动同步；

②心包随着腔室运动；

③忽略心脏瓣膜的运动对心脏腔室的影响；

④将心壁上心肌的方向取为心壁中间心肌的方向，即沿着心壁圆周方向；

⑤心脏腔室和心包表面上的点都沿着法线方向运动；

⑥我们将心动周期分为7个时间段，并合理假定每一时间段内心脏腔室和心包都近似为线性运动。

根据心脏运动的生理特性，我们把心动周期分为七个时间段，分析每个心动周期中各个时间段的腔室压力和容积变化。图9为心动周期各时期压力变化和心室容积变化图，其中C1表示心房收缩，C2表示等容收缩期，C3表示快速射血，C4表示缓慢射血期，C5表示舒张前期，C6表示等容舒张期，C7表示快速充盈相，C8表示缓慢充盈相，C9表示主动脉压力，C10表示心室容积，C11表示心房压力，C12表示心室压力，C13表示心电图。图10是我们将左心室的容积变化线性化的结果。

基于前面的假设和近似，动态心脏建模的大体步骤可以分为三步：首先建立心脏腔室的运动模型，接着建立心包的运动模型，最后将它们组合在一起构成整个心脏的四维动态模型。由于在建立腔室动态模型的过程中已经考虑到了相邻腔室之间的作用，因此不会出现组合在一起发生冲突的情况。最后可将整个心脏的四维动态模型用VTK进行显示：

(1.1.4.1)构建腔室动态模型：腔室动态建模步骤的主要思想就是从腔室体积的变化推出腔室的形变，即$V\left(t\right)\Rightarrow r\left(t\right),$这里V(t)代表腔室的体积随时间的变化，r(t)代表腔室的形变随时间变化(即曲面上的点的位移随时间的变化)。首先根据心脏的静态模型的形态确定它在一个心脏运动周期中所处的位置t₀，然后按照下述步骤分别建立心脏四个腔室的动态模型：

1)根据心脏的运动特性，将心脏一个周期分为7个不同的时期，选取这七个时期的起始和结束的时间点作为7个参考时间点t_i(i＝1，...，7)，其中t₁为t0下一个时间。

2)在时间t_i时，根据曲线V(t)求出相对于时间t_i-1的体积变化ΔV。

3)计算左心室表面(不包括腔室连接处)的曲面面积S。所谓的腔室连接处就是房室瓣和室间隔这些地方。

4)计算从时间t_i-1到t_i中间每个点产生的沿着该点法线方向的位移量为d＝ΔV/S。

5)对曲面上(不包括腔室连接处)的每一个取样点P，计算在点P处的曲面法线方向令$\overrightarrow{{\mathrm{OP}}^{\prime }}=\overrightarrow{\mathrm{OP}}-d·\overrightarrow{n},$其中O为坐标原点，点P′为形变后的取样点的新位置，在左心室中，

$d=\left(\begin{array}{cc}210t/S\left(t\right),& t\in \left[\mathrm{0,0.1}s\right)\\ 0,& t\in \left[\mathrm{0.1,0.15}\right)\\ 470(0.15-t)/S\left(t\right),& t\in \left[\mathrm{0.15,0.25}\right)\\ \frac{460}{3}(0.25-t)/S\left(t\right),& t\in \left[\mathrm{0.25,0.40}\right)\\ 0,& t\in \left[\mathrm{0.40,0.47}\right)\\ \frac{4600}{11}(t-0.47)/S\left(t\right),& t\in \left[\mathrm{0.47,0.58}\right)\\ \frac{150}{11}(t-0.58)/S\left(t\right)& t\in \left[\mathrm{0.58,0.80}\right)\end{array}\right),$

S(t)为在t时刻的心脏腔室表面积。

6)用前面的建模步骤来对新的取样点重构心脏模型。存储新的心脏模型和对应的时间t_i。

7)移到下一个时间点t_i+1，重复上面的第2到第6步直到一个心脏周期完成。

8)将对应于选取时间点t_i的心脏腔室模型连接起来构成腔室的动态模型，令相邻时间点间的心脏腔室上的点近似为线性运动；

(1.1.4.2)心包动态模型构建：心包的建模是基于腔室的动态模型和心肌的收缩模型，步骤在下面列出。在过程中基于另一个假设，即心肌体积不变。令心脏在收缩前、后的体积为V₁和V₂，d₁和d₂是心肌在心包和腔室之间的厚度，L₁和L₂是心肌纤维的长度。根据之前的假设，我们可以得到下列公式：

$\left(\begin{array}{c}{L}_1{d}_1=L_2{d}_2\\ \frac{V_1}{V_2}=\frac{L_1^3}{L_2^3}\end{array}\right)\Rightarrow \frac{d_1}{d_2}={\left(\frac{V_1}{V_2}\right)}^{1/6}$

具体步骤如下：

1)从时间t₀开始。

2)对心包的每一个取样点P，计算在这一点的法线找到法线与腔室的交叉点Q。

3)在时间t＝t_i时，根据腔室动态建模算法得到点Q的位移矢量

4)根据心包和腔室之间的心肌层厚度计算心脏壁的厚度d＝|PQ|，然后得到新的心脏壁厚度$d^{\prime }=d·{\left(\frac{V_i}{V_{i+1}}\right)}^{1/6}.$

5)令$\overrightarrow{{\mathrm{OP}}^{\prime }}=\overrightarrow{\mathrm{OP}}+(\overrightarrow{s}·\overrightarrow{v}+d-d^{\prime })·\overrightarrow{n}$为形变之后的取样点新位置矢量。

6)根据新的取样点重建心包模型，保存模型及其对应的时间t_i。

7)移动到下一个时间点t_i+1，重复建模的步骤2)到步骤6)直到一个心脏周期结束。

8)将对应于选取时间点t_i的心包模型连接起来构成心包的动态模型，令相邻时间点间的心包上的点近似为线性运动；

(1.1.4.3)心脏动态模型的建立：心脏的动态模型是由一个心动周期内一组不同时刻的三维模型构成的。在上面的假设中我们将心脏周期分成了7段，并假设每段里的心脏近似为线性运动。由于我们的原始数据是人体的解剖数据，根据相关资料关于心脏停止跳动后的生理特性以及心脏静态模型中计算出来的体积大小推测前面所建立的原始静态模型处于心房收缩期。因此我们可以先建立8个时间点的瞬态三维模型，然后用线性变化来描述它们之间的过渡过程，这样就构成了心脏整个周期的运动。这八个时间点分别对应着动态模型的初始时刻(即原始静态模型在心脏周期中所处的时刻)，和心脏周期的7个分段的首末时间点。这里需要注意一点的就是相邻时间点之间的时间间隔并不是相等的，而是根据实际心动周期的时相来确定的。设置一个定时器，来响应不同时刻的心脏模型数据，可以把动态心脏在VTK中显示出来。

图11是我们将建立的周期包含7个阶段的心脏动态模型用VTK显示的结果，每幅子图各为每个阶段的一个实例。

(1.2)冠状动脉四维动态模型建立

建立冠状动脉四维动态模型是为了指导从多幅不同投影角度X射线造影图中的冠状动脉血管三维重建，其步骤是：

(1.2.1)血管切片轮廓提取步骤：根据解剖学的知识和人体横断解剖的图片，在每幅原始图片中分出血管。图12是原始心脏切片图像及其血管提取结果对比图。

(1.2.2)提取血管中轴点与半径步骤：建立冠状动脉模型要提取血管的中轴点和半径两个信息。中轴点取定为血管切片轮廓的形心，同时假设血管为圆形，血管的半径为轮廓的内切圆半径。

(1.2.3)重建血管三维骨架步骤：把上一步提取的冠状动脉各个分支的中轴点用B样条曲线进行拟合得到血管的三维骨架模型。

(1.2.4)带半径的血管三维模型的构建步骤：在所述的血管三维骨架模型上使用GC模型加入半径信息，建立静态的冠状动脉模型。

广义圆柱体表示又叫广义锥表示，它是一种推扫模型。由轴线、截面和扫描规则组成。轴线是一条三维空间曲线，截面是一个平面区域，它沿着轴线以一定角度进行扫描，截面的形状在扫描过程中可以变化，从而可以表示不同的形体。扫描规则定义了截面与轴线之间的几何关系。物体的表面就是这些横截面封闭曲线的并集，推扫形成的体积就是广义圆柱体。这种表示能很好的描述许多具有轴对称性的物体，相似的物体有相似的轴和横截面。

广义圆柱(GC)用一个三元式(A，E，α)来定义血管。A(X，Y，Z)(s)是轴线，它是由参数s定义的一条空间曲线；α是倾斜角，它是由横截面所在平面和在A(s)处轴线切线的夹角；E＝(t，s)是一条平面曲线，它以在A(s)的平面的参数定义了横截面，可以分为两个函数：E(s，t)＝r(s)C(t)，其中轮廓函数C(t)描述了横截面的形状，而半径函数r(s)描述了它的尺寸。这里我们假设血管的横截面是圆形C(t)＝(cos2πt，sin2πt)，并垂直于中轴的切线方向(即倾斜角α为90度)。

利用此模型对血管进行三维体描述，认为血管树是一系列横截面的形状和大小以及中轴线的方向都不断变化的管状(圆柱)结构的集合。

(1.2.5)最后对模型中的各个血管基元按照Dodge et al.(Dodge JT，Brown BG，Bolson EL，Dodge HT：“Intrathoracic spatial location of specified coronary segments on thenormal human heart.”Circulation 78：1167-1180，1988.)设计的命名法进行命名。图13是我们对模型中血管基元命名的结果。

(1.2.6)建立冠脉四维动态模型步骤：冠状动脉的运动是由心脏的运动引起的，它与心脏一样进行周期性运动。假设血管保持管状，在运动过程中半径不发生变化。将血管的运动建模分成两步。

1)建立血管的中轴的运动模型

采样血管中轴上面的一些点，使它们按照心包运动的规律运动，根据心脏生理学知识，心动周期分为7个阶段，合理地假设这些取样点在这7个阶段里分别近似为线性运动。对血管取样点通过分别建立这7个阶段内的运动模型，得到它们在整个心脏周期的的运动模型。然后对这些取样点在同一时刻的位置进行B样条曲线拟合得到对应于每一时刻的血管中轴模型，一个周期中所有时刻的血管中轴模型就组成了血管中轴的周期运动模型。图14是建立的一组七个时刻的血管中轴运动模型。

2)在中轴的运动模型上加半径信息

这里同样通过GC模型在动态的血管骨架模型上加入半径信息，得到动态的冠状动脉模型。

(2)从不同角度的至少二幅X-射线造影图进行血管三维重建(如图3所示为血管三维重建流程图)：

(2.1)血管造影图预处理：为了消除干扰噪声，保留并突出所需的图像信息，最终准确地从造影图像中提取血管结构，首先对图像序列采用频域增强的方法来增强血管，然后采用形态学Bottom-Hat变换来进一步强调血管。

(2.1.1)频域增强

由于心脏是有运动周期的，由造影系统成像得到的为一个冠状动脉造影图像序列。研究这一图像序列，我们发现常见的肋骨、脊髓、肺和一些内部组织可以看作是不变或者缓慢变化的信号，然而冠状动脉随着心脏的收缩和扩张，呈现与心脏有着相似频率的运动状态。图15是左冠状动脉血管和背景噪声点的灰度在造影图序列中的变化趋势，其中选取的是北京朝阳医院的一个病人的造影图序列。其中曲线1表示坐标为(252，134)的血管像素的灰度值，曲线2表示坐标为(362，95)的背景点的灰度值。从图中可以看出血管点的灰度变化比背景噪声点要大，根据这个特点可以通过高通滤波去掉背景噪声。

设一帧数为T的图像序列f₀(x，y，t)，单幅图像大小为M×N，0≤x＜M，0≤y＜N，0≤t＜T，x，y，t均为整数。对图像序列采用如下方法进行频域增强：

1)对图像序列f₀(x，y，t)进行M×N次沿时间轴t的离散傅立叶变换，得到频域信号F₀(x，y，k)，0≤k＜T，离散傅立叶变换公式如下所示：其中j为虚数单位

$F_0(x,y,k)=\frac{1}{T}\underset{t=0}{\overset{T-1}{\Sigma }}{f}_0(x,y,t)e^{-j2\mathrm{\pi kt}/T},$

0≤k＜T，0≤x＜M，0≤y＜N

2)对频域信号F₀(x，y，k)采用下面的式子进行M×N次沿k轴的高通滤波，

F(x，y，k)＝F₀(x，y，k)(1-e^-βm)，

高通滤波函数为H(x，y，k)＝(1-e^-βt)，其中β是预先给定的常数，用以抑制某些频率的背景噪声。理论上，β的取值范围是[0，∞]，当β→0时，(1-e^-βt)→0，F(x，y，k)→0，所有的信号包括冠状动脉都被去除，使得抑制噪声过度；当β→0，(1-e^-βt)→1，F(x，y，k)→F₀(x，y，k)所有的信号包括冠状动脉都被保留，没有起到去除噪声的作用。合适的β能够最大程度上去除噪声并保持冠脉信号。经过多次实验，β值的优选范围为1/8～1/5。

3)利用下式离散傅立叶反变换作用于滤波后的频域信号，获得新的图像序列：

$f(x,y,t)=\underset{k=0}{\overset{T-1}{\Sigma }}\bar{F}(x,y,k)e^{j2\mathrm{\pi k}/T}$

(2.1.2)多尺度形态学Bottom-Hat变换

我们采用形态学Bottom-Hat变换来进一步增强血管。Bottom-Hat变换是一种常见的用于提取图像暗结构的形态学算子，定义为g＝f-(f·b)，其中，f是输入图像，b是结构元素函数，$f·b=(f\oplus b)\mathrm{\Theta b}$代表闭操作，由膨胀操作和腐蚀Θ操作组成。闭操作可以去除由结构元素b表示的图像中较暗的感兴趣结构，并保留较亮像素。在造影图像中，目标物体(血管)的灰度值较背景低，因此，f·b的结果可以看作是背景，由原图减去背景图像可以得到血管部位。

在实现中，我们选择两个不同直径尺度D₁，D₂(D₁＞D₂)的圆盘结构算子，较大直径D₁需略大于最宽的冠状动脉，较小直径D₂略大于最窄冠状动脉以集中选取对比度不高的小动脉。我们对增强图像f(x，y)分别进行尺度为D₁和D₂的Bottom-Hat变换，得到g₁(x，y)和g₂(x，y)。

经过上述变换后得到的图象是突出血管的灰度级图像，但是仍然含有背景噪声，需要分别对g₁(x，y)和g₂(x，y)进行以下步骤以得到二值图像。

1)相对阈值比较。Bottom-Hat变换强化了图像中的暗感兴趣结构并减弱了其他部分，使得Bottom-Hat变换后造影图的血管部分像素值高于变换前造影图中的血管部分像素值，背景部分像素值低于变换前的背景部分像素值。通过选取合适的阈值T₁，将Bottom-Hat变换之前图的像素值和Bottom-Hat变换之后图的像素值与T₁进行比较，可以得出哪些点是目标，哪些点是背景。由于图像之间的差异很大，根据图像灰度值分布自动地计算阈值有着更广泛的适用性。我们把这种方法称为“相对阈值比较”，它能够得到将血管与背景分割开的二值图：

首先，采用下式比较g(x，y)和f(x，y)得到差值图d(x，y)，

$d(x,y)=\left(\begin{array}{cc}g(x,y)-f(x,y)& \mathrm{if}(g(x,y)-f(x,y)>0)\\ 0& \mathrm{otherwise}\end{array}\right),$

然后，根据d(x，y)自适应计算阈值T₁，T₁为d(x，y)图中的非零像素值的平均值，

$T_1=\frac{\underset{d(x,y)>0}{\Sigma }d(x,y)}{N_d},$

其中N_d是d(x，y)中非零像素值的个数。

最后，若d(x，y)≥T₁，则该点为血管象素点，若d(x，y)＜T₁，则该点为背景点。将g₁(x，y)和g₂(x，y)分别进行上述处理得到二值图B₁(x，y)和B₂(x，y)。

2)整合二值图。叠加二值图B₁(x，y)和B₂(x，y)，目的是加强末端血管和小血管，同时去除一些面积低于阈值T₂的块状噪声区域。此时，我们就得到了冠状动脉整个提取结果。

图16是将原始血管造影图经过频域增强和Bottom-Hat变换之后的结果。图16(a)是原始造影序列图的第22帧，图16(b)是经过频域增强之后的结果，图16(c)是将增强图像进行Bottom-Hat变化之后的结果。

(2.2)提取血管骨架及半径：对预处理后的图像进行分割提取血管区域，采用形态学的方法细化血管区域得到血管骨架。选用八连通链码来进行心血管骨架的跟踪，提取血管半径，用二叉树存储血管树的拓扑结构。

(2.2.1)通过八连通链码来进行心血管骨架的跟踪：

Step1.搜索起始点从边缘图的左上角开始，由上往下，由左至右方向搜索。取dir＝7，并保存一条新轮廓的起始坐标(r，c)。(其中dir为扫描方向变量，记录上一步中沿着前一个边界点到当前边界点的移动方向，下表列出了不同的dir值对应的不同的移动方向)。

  dir  0  1  2  3  4  5  6  7  X向变化  1  1  0  -1  -1  -1  0  1  Y向变化  0  1  1  1  0  -1  -1  -1

Step2.按逆时针方向搜索当前像素的3×3邻域，其起始搜索方向设定如下：

若dir为奇数，取(dir+7)mod 8；

若dir为偶数，取(dir+6)mod 8。

在3×3邻域中搜索到的第一个与当前像素值相同的像素便为新的边界点A_n，同时更新变量dir为新的方向值。在对当前边缘点进行链码编码后，将当前边缘点的像素置为背景像素值，使下次检测边缘不会再搜索到。

Step3.如果检测到一个边缘点的3×3邻域已经没有边缘点存在，则该点是轮廓的末端点。随后储存在相应的轮廓链码中。

在提取血管骨架之后，根据下列步骤测量心血管直径：用圆盘模板沿着骨架方向寻找直径，在需要提取直径的骨架点处，圆盘模板的直径在适当的范围内增大，当圆盘与边缘点只要出现一个相切点时，此时圆盘模板的直径为该点血管直径。利用血管直径变化的连续性，借助骨架链码，沿着骨架的方向在一定范围内增大或缩小血管直径，可减少搜索范围，提高搜索效率。

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

用树型结构来表示血管的拓扑结构，首先仅需手动指定血管的根结点的位置，然后通过二叉树搜索得到血管的树型结构表示。在这个过程中会出现两个问题：第一，由于血管造影图是通过X射线透射人体得到的，通过投影三维空间中的不同血管可能会在投影图上形成交叉，必须要能够识别出其中的血管交叉，否则二叉树搜索得到的结果并不是正确的血管基元投影；第二，由于实际X射线造影过程中的各种因素影响以及人体组织对X射线产生的衰减，会在造影图中引入很多噪声，从而可能导致伪血管的出现，这个也是必须要消除的。

(2.2.2.1)血管交叉的识别

理想情况下血管交叉于一点，如果血管中轴上点连接的血管段数为4，就判断血管在这一点出现交叉。若只出现这种情况，也可以很容易地去除交叉，即先将血管段断开，再把相对应的血管段连接起来。然而实际情况下，在血管造影图的血管中轴提取图像中，血管的交叉部分可能是一小段血管而不是一个点，这样就不能简单地用上面的方法来判断了。

从血管造影图的血管中轴提取图像中可以看到，当出现血管交叉时，二叉树搜索会出现循环。图17是血管交叉时二叉树搜索出现循环的示意图。

利用这一特点可判断，当二叉树搜索过程中出现循环时，则说明有血管发生了交叉。

在发现搜索过程中出现循环后，下面需要找到发生交叉的血管段：

记录发生循环的血管段组成的环，根据二叉树搜索的特点，交叉血管(即两个血管段的交叉部分)应该在环上或连接在环上。逐个搜索环上的血管或连接在环上的血管结构(即连接在环上的血管段上但不在环上)，当它们满足以下条件时，将其作为交叉血管的候选：

Number(children)≥2，

即当它的子血管数大于或等于2才有可能是交叉的血管基元投影，设置代价函数

f＝a₁·Length(seg)+a₂·min{fabs(angle(seg1，seg3)-π)+fabs(angle(seg2，seg4)-π)，

                           fabs(angle(seg1，seg4)-π)+fabs(angle(seg2，seg3)-π)}

其中a₁，a₂为加权系数，根据血管基元长度和血管分支角度的归一化系数选取。seg为当前的候选血管结构，seg1和seg2为环上连接seg的血管结构，seg3和seg4为seg另一边连接的血管结构，如图18所示。

选取其中代价函数最小的候选血管结构作为交叉血管基元投影。然后根据下面的函数进行平滑性判断，

g₁＝fabs(angle(seg1，seg3)-π)+fabs(angle(seg2，seg4)-π)

g₂＝fabs(angle(seg1，seg4)-π)+fabs(angle(seg2，seg3)-π)

若g₁＜g₂，则seg1和seg3属于同一条血管，seg2和seg4属于同一条血管；若g₁≥g₂，则seg1和seg4属于同一条血管，seg2和seg3属于同一条血管。

将属于同一条血管的两段血管连接起来得到一条完整的血管，从而去除投影形成的血管交叉。

(2.2.2.2)伪血管的识别与去除

本发明中讨论的伪血管主要针对在造影图中出现的，由于人体本身组织结构(如肋骨)、造影剂分布不均匀或者注入造影剂过程中的导管和噪声等原因，而在图像中留下了一些类似血管而非血管的图形。分割过程中将这些图形误认为血管分割出来，形成了本章节中提到的伪血管。

由于本发明中的模型只是表达血管的基本结构，没有也不大可能考虑到个体的差异，因此无法用模型来判别伪血管。于是利用血管的拓扑结构关系来判别伪血管：一般情况下，伪血管出现在两个角度的造影视图中时拓扑结构是不相同的；此外根据投影几何知识，投影并不会改变血管的拓扑结构(这里是指定性而非定量的连接关系)。基于上述条件，通过比较两幅造影图中是否存在拓扑结构不同的血管结构来判断此血管结构是否为伪血管，如图19所示。

(2.2.2.3)实现步骤

考虑到上面的问题，本发明设计了如下步骤：

1)根据二叉树搜索的方法从血管根结点开始搜索，并对搜索过的血管像素做标记，当搜索到标记过的血管像素时则停止向前搜索，这样就可避免发生死循环；

2)搜索整个树，检测是否出现循环，若出现循环，则按照3.1.1中提到的方法找出交叉血管结构并纠正；

3)重复步骤2，直到没有循环出现为止。

(2.3)血管造影图标记

通过匹配血管造影图与血管树模型投影图中的血管，对血管造影图中的血管进行标记。这里的血管树模型投影图是将血管三维模型仿照血管造影系统投影的方式，采用形成相应血管造影图中的参数进行投影得到的，因此可以通过比较两段血管之间的拓扑相似性来判断它们是否匹配。

由于个体差异性，血管造影图与血管树模型投影中的血管数量和血管分布并不相同。由于建立的血管树模型是只取了冠脉系统比较显著的血管，能够代表大多数人的共性，因此合理地假设血管树模型里的血管在造影图中都有对应的血管。

首先需要选定用来指导造影图中血管标记的血管树模型的尺度，由于所选病人造影图中分割得到的血管并不多，这里选用二级血管树模型对它进行标记。通过匹配血管造影图和血管树模型中的血管结构，对血管造影图中的血管进行标记。

血管标记步骤：

(I)按X射线造影方式从不同方向对血管树模型进行投影得到不同视角的投影图。

(II)求每个投影图与血管造影图的血管中轴提取图像T的相似度，找到其中相似度最高的投影图P。

相似度通过以下步骤求取：

a、将血管树模型投影图和血管造影图的血管中轴提取图像T中的血管拓扑结构分别保存在二叉树T1和T2中；

b、从T2中提取与T1拓扑结构匹配的子树T3。对T1和T3进行匹配，计算下列代价函数：

f＝a₁·f₁(angle)+a₂·f₂(curve)+a₃·f₃(length)

其中f₁(angle)表示两者之间血管拓扑结构的相似性，表示为血管分支之间角度angle的函数；

f₂(curve)表示两者之间血管结构曲线形状的相似性，表示为曲线曲率curve的函数；

f₃(length)表示两段血管长度的相似性，表示为血管长度length的函数；

a₁，a₂，a₃为这三种相似性的加权系数，根据血管分支的角度、血管曲率和血管长度的归一化系数选取。取其中代价函数最小的血管对为匹配血管。

c、取其中代价函数最小的子树为T1的匹配子树，将此最小代价函数的倒数作为模型投影图与血管造影图之间的相似度。

(III)通过上述求相似度的步骤得到与血管造影图的血管中轴提取图像T相似度最高的投影图P在T中的匹配子树M。用P中血管基元投影的名字对子树M中的血管基元投影进行标记。图20是进行标记的初步结果，(a)图中的G1，G2，...，G7分别与(b)图中的H1，H2，...H7标记相同。

(IV)当选取的两个造影图不在心动周期中的同一时刻时，要分别用这两个时刻的血管树模型对相应造影图的血管基元投影进行标记。

(2.4)血管匹配

血管匹配过程包括以下两个步骤：

1)找到在两幅造影图中有相同标记的一对血管基元投影作为匹配的血管基元投影；

2)对照两幅造影图，根据拓扑结构的限制在两幅造影图中没有标记的血管基元投影里找到拓扑结构相同的血管基元投影作为匹配的一对血管基元投影；

对于两幅造影图中还没有标记的血管结构，可以通过以下方法找匹配的血管对：

a、设这两幅造影图中血管拓扑结构分别存在left_T和right_T中。

b、找到left_T中没有标记的血管UnName1，在right_T中相同位置处搜索是否也有未标记的血管结构UnName2，如果有就说明可能存在对应的血管结构，没有则判断UnName1为伪血管；

c、若UnName2存在，需判断UnName1和UnName2之间是否一一对应。考虑到有些情况下，造影图中相同位置处可能存在不止一个未标记血管结构，又要分为两种情况进行讨论：若是一一对应，则直接将UnName1和UnName2命以相同标记；若非一一对应，如图21中所示，则需要两两之间进行匹配，找到代价函数最小的血管对为匹配血管，并作相同标记。

d、继续步骤a，b，c，直到血管树里的血管全部标记或判断为伪血管为止。

e、图22是对两幅造影图最后的匹配结果，其中造影图22(a)中的血管段J1，J2，...，J13分别与22(b)中的血管段K1，K2，...，K13匹配对应。

通过对很多个病人的造影图进行实验，三级血管以内的大血管都可以正确标记，四级以下小血管标记准确率不是很高，这与图像质量和模型的细节都有关系。随着模型细节的进一步补充，以及模型的进一步优化，小血管标记的准确率将会逐渐提高。

(5)三维重建：进一步寻找匹配血管基元投影上各像素点之间的对应，根据重建公式对血管进行三维重建。

(5.1)在重建之前，有一点要先考虑：在进行血管三维重建时，选取的两幅造影图必须对应同一时刻。而我们得到的图像一般是通过单臂造影系统得到的一个造影序列图像，在单臂造影系统中的造影图像都对应不同的时间。为了满足选取的两幅造影图在统一时刻的要求，研究者们大都采用选取舒张末期的血管造影图来做三维重建，但是这样会浪费掉其他时刻的造影图片，可能会丢失很多信息。而且这样也没有考虑到造影图的效果，比如血管遮挡、噪声等，很难保证选定时刻的这些不同方向的血管造影图适合于做三维重建。如果能够使其他时刻的血管造影图都可用，将会大大增加可用信息量，提高重建的精度。

当选取的两幅造影图不在同一时刻时，由于它们并不对应一个静止的心脏，不能用一般的三维重建公式，否则会有大的误差。因此需要先将这两幅图中的血管尽可能统一到同一个时刻，即进行心脏运动补偿：

假设现在选取了两幅造影图I₁和I₂，分别对应于心动周期中的t₁，t₂时刻，I₂的投影角为θ。设I₂中的血管像素为P，t₁时刻在θ方向的造影图为I₃，P在I₃中的对应点为P’，由P变换到P’的过程就是运动补偿。

设心脏模型在t₁，t₂时刻在角度θ上的造影图分别为J₁，J₂。点P在J₁中的对应点为Q，点P’在J₂中的对应点为Q’，O为原点。则可由以下公式计算运动补偿：

$\overrightarrow{{\mathrm{OP}}^{\prime }}-\overrightarrow{\mathrm{OP}}=\mathrm{scale}·(\overrightarrow{{\mathrm{OQ}}^{\prime }}-\overrightarrow{\mathrm{OQ}}),$

其中scale为将造影图中血管归一化到血管树模型时的归一化系数。

将I₂通过运动补偿后，得到血管新的位置I₃。利用I₁和I₃进行血管三维重建，其中I₁和I₂的血管匹配关系对I₁和I₃仍然适用。

(5.2)基于对应点的三维重建

本节利用两个不同角度的单面冠脉造影图像上的离散点重建三维冠状动脉三维离散点，简单地说，也就是根据两幅不同角度的投影图像上的对应点，计算出空间点的三维坐标。

(5.2.1)冠脉造影系统三维空间上的坐标变换

附图27是造影系统在两个不同角度的成像示意图。利用坐标变换的基础知识，图27中从坐标系x₁y₁z₁s₁到x₂y₂z₂s₂的几何变换运动是三维空间的刚体运动，可以用旋转运动和平移运动描述。而且运动过程中x₁y₁s₁平面与s₁z₁轴的相对位置没有变化，所以只需通过坐标变换，使s₁z₁轴和s₂z₂轴重合，即可使坐标系x₁y₁z₁s₁和x2y₂z₂s₂重合。

为方便描述，做以下约定：

(1)(x_1i，y_1i，z_1i)和(x_2i，y_2i，z_2i)分别表示点P_i在三维坐标系x₁y₁z₁s₁和x₂y₂z₂s₂中的坐标；

(2)(u_1i，v_1i)和(u_2i，v_2i)分别表示点P_i在投影平面坐标系U₁V₁O₁和U₂V₂O₂中的坐标；

(3)X射线源S₁和S₂与旋转中心(isocenter)之间的距离分别为L₁和L₂。

由XYZO坐标系转换到x₁y₁z₁s₁坐标系，使ZO轴和s₁z₁重合，可以通过两次旋转变换、一次平移变换实现：(1)绕OY轴顺时针旋转α₁；(2)绕OX轴逆时针旋转β₁；(3)从O点平移到s₁点。

即

[x₁，y₁，z₁]^T＝R_X(β₁)·R_Y(-α₁)·[X，Y，Z]^T+T₁                    (1)

所以

${[X,Y,Z]}^T=R_Y^{-1}(-{\alpha }_1)·R_X^{-1}\left({\beta }_1\right)·({[x_1,y_1,z_1]}^T-T_1)---\left(2\right)$

同样可以通过两次旋转变换和一次平移变换实现XYZO坐标系到x₂y₂z₂s₂坐标系的转换：(1)绕OY轴顺时针旋转α₂；(2)绕OX轴顺时针旋转β₂；(3)从O点平移到s₂点。即

[x₂，y₂，z₂]^T＝R_X(β₂)·R_Y(-α₂)·[X，Y，Z]^T+T₂                    (3)

所以：

${[X,Y,Z]}^T=R_Y^{-1}(-{\alpha }_2)·R_X^{-1}\left({\beta }_2\right)·({[x_2,y_2,z_2]}^T-T_2)---\left(4\right)$

综合式(2)和(3)得：

[x₂，y₂，z₂]^T＝R([x₁，y₁，z₁]^T-t)                                  (5)

其中$R=R_X\left({\beta }_2\right)R_Y\left({\alpha }_2\right)R_Y^{-1}\left({\alpha }_1\right)R_X^{-1}(-{\beta }_1)t=-R^{-1}{T}_2+T_1$

T₁＝(0，0，L₁)^T   T₂＝(0，0，L₂)^T

(5.2.2)点的3D重建方法

在图27中，假设坐标系x₁y₁z₁s₁中冠状动脉树骨架点P的三维坐标设为(x_1i，y_1i，z_1i)，P在图像A坐标系U₁V₁O₁上的投影点p₁坐标为(u_1i，v_1i)；坐标系x₂y₂z₂s₂中P的三维坐标设为(x_2i，y_2i，z_2i)，P在图像B坐标系U₂V₂O₂上的投影点p₂坐标为(u_2i，v_2i)。D₁和D₂分别表示X射线源S₁和S₂到各自投影平面的垂直距离。

若已知(u_1i，v_1i)和(u_2i，v_2i)，由下述推导，可以根据(R，t)求解P的三维坐标(x_1i，y_1i，z_1i)：

根据透视投影的几何关系，可知

$\left(\begin{array}{c}{x}_{1i}\\ y_{1i}\\ z_{1i}\end{array}\right)=z_{1i}·\left(\begin{array}{c}{\xi }_{1i}\\ {\eta }_{1i}\\ 1\end{array}\right),$$\left(\begin{array}{c}{x}_{2i}\\ y_{2i}\\ z_{2i}\end{array}\right)=z_{2i}·\left(\begin{array}{c}{\xi }_{2i}\\ {\eta }_{2i}\\ 1\end{array}\right)---\left(6\right)$

其中

${\xi }_{1i}=\frac{u_{1i}}{D_1},$${\eta }_{1i}=\frac{v_{1i}}{D_1}---\left(7\right)$

${\xi }_{2i}=\frac{u_{2i}}{D_2},$${\eta }_{2i}=\frac{v_{2i}}{D_2}---\left(8\right)$

进一步可以得到

$\frac{x_{1i}}{z_{1i}}=\frac{u_{1i}}{D_1}={\xi }_{1i},$$\frac{y_{1i}}{z_{1i}}=\frac{v_{1i}}{D_1}={\eta }_{1i},$$\frac{x_{2i}}{z_{2i}}=\frac{u_{2i}}{D_2}={\xi }_{2i},$$\frac{y_{2i}}{z_{2i}}=\frac{v_{2i}}{D_1}={\eta }_{2i}---\left(9\right)$

由式(2)可以得到

$\left(\begin{array}{c}{x}_{2i}\\ y_{2i}\\ z_{2i}\end{array}\right)=R(\left(\begin{array}{c}{x}_{1i}\\ y_{1i}\\ z_{1i}\end{array}\right)-\overrightarrow{t})=\left(\begin{array}{c}(r_{11}·x_{1i}+r_{12}·y_{1i}+r_{13}·z_{1i})-(r_{11}·t_1+r_{12}·t_2+r_{13}·t_3)\\ (r_{21}·x_{1i}+r_{22}·y_{1i}+r_{23}·z_{1i})-(r_{21}·t_1+r_{22}·t_2+r_{23}·t_3)\\ (r_{31}·x_{1i}+r_{32}·y_{1i}+r_{33}·z_{1i})-(r_{31}·t_1+r_{32}·t_2+r_{33}·t_3)\end{array}\right)---\left(10\right)$

$R=\left(\begin{array}{ccc}{r}_{11}& r_{12}& r_{13}\\ r_{21}& r_{22}& r_{23}\\ r_{31}& r_{32}& r_{33}\end{array}\right)$$\overrightarrow{t}=\left(\begin{array}{c}{t}_1\\ t_2\\ t_3\end{array}\right)$

联立式(9)(10)可以解得：

$\left(\begin{array}{ccc}1& 0& -{\xi }_{1i}\\ 0& 1& -{\eta }_{1i}\\ r_{11}-r_{31}·{\xi }_{2i}& r_{12}-r_{32}·{\xi }_{2i}& r_{13}-r_{33}·{\xi }_{2i}\\ r_{21}-r_{31}·{\eta }_{2i}& r_{22}-r_{32}·{\eta }_{2i}& r_{23}-r_{33}·{\eta }_{2i}\end{array}\right)·\left(\begin{array}{c}{x}_{1i}\\ y_{1i}\\ z_{1i}\end{array}\right)=\left(\begin{array}{c}0\\ 0\\ \overrightarrow{a}·\overrightarrow{t}\\ \overrightarrow{b}·\overrightarrow{t}\end{array}\right)---\left(11\right)$

其中：

$\overrightarrow{a}=\left(\begin{array}{ccc}{r}_{11}-r_{31}·{\xi }_{2i}& r_{12}-r_{32}·{\xi }_{2i}& r_{13}-r_{33}·{\xi }_{2i}\end{array}\right)---\left(12\right)$

$\overrightarrow{b}=\left(\begin{array}{ccc}{r}_{21}-r_{31}·{\eta }_{2i}& r_{22}-r_{32}·{\eta }_{2i}& r_{23}-r_{33}·{\eta }_{2i}\end{array}\right)---\left(13\right)$

式(10)可以简写为

A·C＝B

$A=\left(\begin{array}{ccc}1& 0& -{\xi }_{1i}\\ 0& 1& -{\eta }_{1i}\\ r_{11}-r_{31}·{\xi }_{2i}& r_{12}-r_{32}·{\xi }_{2i}& r_{13}-r_{33}·{\xi }_{2i}\\ r_{21}-r_{31}·{\eta }_{2i}& r_{22}-r_{32}·{\eta }_{2i}& r_{23}-r_{33}·{\eta }_{2i}\end{array}\right)---\left(14\right)$

$C=\left(\begin{array}{c}{x}_{1i}\\ y_{1i}\\ z_{1i}\end{array}\right)$$B=\left(\begin{array}{c}0\\ 0\\ \overrightarrow{a}·\overrightarrow{b}\\ \overrightarrow{b}·\overrightarrow{t}\end{array}\right)$

由式(14)可以看出，该方程组由4个线性方程组成，求解3个未知量x_1i、y_1i、z_1i，因此是一个超限定方程组，其最小二乘解为：

C＝(A^T·A)^-1·A^T·B                               (15)

即坐标系x₁y₁z₁s₁中P点的三维坐标(x_1i，y_1i，z_1i)。根据式(15)，就可以重建出三维坐标点了。

考虑到外极线约束的缺陷，利用提取特征点的方法，对两幅造影图中匹配的一对血管基元投影上的各点进行匹配，得到各匹配像素点对，并进行三维重建得到各空间三维点。用B样条拟合这些重建得到的空间三维点得到三维血管树。

附图23中(g)和(h)分别显示了用模型指导重建算法和外极线匹配算法重建得到的血管造影图三维重建结果，(i)是在两幅造影图中手动取对应点得到的重建结果。

附图24中(a)和(b)是将用本文方法重建的血管三维结构即图23(g)分别按图23(a)和(b)的造影角度进行反投影的结果，(c)和(d)是对用来进行重建的两幅原始造影图提取的中轴图像，分别与图23(c)和(d)相同。通过对比图(a)和(c)，(b)和(d)可以看到重建的结果还是比较好的。

(2.6)更新模型：在重建出血管三维结构后，用得到的正确三维重建结果对模型进行完善。包括以下两个方面的完善措施：

(2.6.1)将三维重建结果与模型进行融合，使模型更加的通用化。这里融合包括以下几个步骤：

1)在重建得到的三维血管结构中找到与模型中血管对应的血管基元为参考血管基元。由于在用模型指导造影图中血管基元的投影标记时得到过血管基元的投影与模型中血管的投影的对应关系，模型上的血管基元在重建得到的三维血管结构上对应的血管基元也因此可以确定；

2)以模型为参考，对参考血管基元进行仿射变换与模型中对应血管基元进行匹配；

3)将仿射变换后的参考血管基元与模型中的对应血管基元上的对应点的坐标取平均得到新的血管树模型。

(2.6.2)对模型细节进行补充。当出现模型中不存在的血管结构时，如果它在对不同病人的造影图进行重建的结果中都出现了，则将它添加到模型中，并将模型中所有的血管按级数来存储，比如一级血管，二级血管等等。由此建立一个血管树的多尺度模型。根据不同的情况，选择不同的尺度模型。

附图25是我们建立的左冠状动脉的多尺度模型示意图。其中：图25(a)是一级血管树模型示意图；图25(b)是二级血管树模型示意图；图25(c)是三级血管树模型示意图；

附图26是我们用重建结果对模型进行补充细节的示意图。其中，(a)是我们由附图23(a)(b)的两幅造影图重建得到左冠状动脉血管三维结构，与图23(g)不同的是图26(a)使用了三级血管树模型对血管进行指导重建。图26(b)是对图26(a)中血管中轴加上半径的结果，(c)是根据虚拟人切片数据建立的血管树模型，(d)是用(a)中重建结果对(c)中模型进行补充细节之后的结果，其中圆圈标出的血管段为补充的血管基元。