CN101901498B - 一种人体冠状动脉血管的递推建模方法 - Google Patents

Abstract

人体冠状动脉血管的递推建模方法，属于数字图像处理与医学成像的交叉领域。本发明首先在由心脏切片图像建立的心脏和冠状动脉血管静态模型的基础上，同时考虑心脏的法向和切向的运动以建立它们的运动模型，然后进一步将从实际单臂X射线造影图像数据重建得到的个体三维冠状动脉血管树与冠状动脉血管模型作差异比较，若拓扑差异较大，则构建新的特异性冠状动脉血管模型，反之，则对两者进行融合，以获得更具一般性的冠状动脉血管模型。本发明可以通过递推的方式不断更新模型以得到一个更加接近人体生理实际的冠状动脉血管模型以及其他具有特异性的冠状动脉血管模型，为临床治疗或研究提供客观的指导和参考。

Description

一种人体冠状动脉血管的递推建模方法

技术领域

本发明属于数字图像处理与医学成像的交叉领域，具体涉及一种人体冠状动脉血管的递推建模方法。

背景技术

国内外在人体心血管建模这方面做的不是很深入，最常用到的一个冠状动脉血管模型是Dodge依据很多个病人的造影图像重建得到的冠状动脉三维结构，建立一个具有共性的冠状动脉血管模型。但是它是一个静态的模型，而且只考虑到了不同个体冠状动脉的共性，不能表达个体的特异性。至于在动态冠状动脉建模方面，也并未明确强调模型的建立，其主要是通过从医学图像(如X射线造影图)序列中提取出血管的运动信息，再将该信息作为可变模型的约束条件，并作用于被初始化的三维静态冠状动脉树使之发生形变，在一定程度下可以理解为动态血管模型。但由于它受到所选取医学图像的限制，因而仅能反映出个体的特异性，而不能代表共性。显然，如何建立既具有共性又具有个性的冠状动脉血管模型是一个难题。

另外，国外也有人做过心脏动态模型的研究，且以双心室建模较为常见。一般做法是，根据心脏的电生理学建立心脏的电生理方程，利用心脏的心肌模型，通过电传导建立心脏的运动。但是由于很难直接得到人体的心脏的心肌模型，大部分这种心脏动态模型都是通过一些动物(如狗等)的数据建立的，不能表达人体冠状动脉血管的特殊性，更不能满足指导X射线血管造影三维重建的医学成像诊断需求。

申请人于2008年5月28日提交了“动态模型指导下的血管造影三维重建方法”中国专利文献(公开号为CN101301207A，公开日为2008年11月12日)。该方法通过对提取的每层心脏切片轮廓采样，获取心脏表面的三维采样点，再利用B样条对这些采样点进行曲面拟合，建立心脏三维静态模型。而之后在建立心脏动态模型时只考虑到腔室表面法向的运动，未考虑切向运动，这便导致了在心脏动态模型基础上建立起来的冠状动脉血管三维动态模型也存在着同样运动信息不准的缺陷，进而影响到血管造影三维重建的效果。此外，在冠状动脉血管模型更新部分的差异性判据仅利用到个体与个体之间对应血管间的平均距离，而未结合个体在生理坐标系下的分布情况，所以，在有些情况下，并不能准确的反映出实际差异。

发明内容

本发明提供了一种人体冠状动脉血管的递推建模方法，利用该方法可以得到一个更加接近人体生理实际的动态的冠状动脉血管模型。

本发明提供的人体冠状动脉血管的递推建模方法，该方法包括下述步骤：

第1步建立冠状动脉血管三维静态模型，其过程为：

(1.1)提取冠状动脉血管切片轮廓；

(1.2)提取冠状动脉血管的中轴点与半径，中轴点为冠状动脉血管切片轮廓的形心，同时假设冠状动脉为圆形，冠状动脉血管的半径为轮廓的内切圆半径；

(1.3)重建冠状动脉血管三维骨架模型：把步骤(1.2)提取的冠状动脉血管各个分支的中轴点用B样条曲线进行拟合得到冠状动脉血管三维骨架模型，即初始时刻下的冠状动脉血管中轴三维静态模型；

(1.4)构建带半径的冠状动脉三维模型：在所述的冠状动脉血管三维骨架模型上使用广义圆柱模型加入半径信息，建立带半径的冠状动脉三维模型，即冠状动脉血管三维静态模型；

第2步建立冠状动脉血管三维动态模型，其过程为：

(2.1)建立冠状动脉血管中轴三维动态模型：

利用冠状动脉血管与心包的附着关系，将心包动态模型的运动信息转化为冠状动脉血管的运动信息，该运动信息包括法向运动矢量和切向运动矢量；在冠状动脉血管中轴取等间隔的点作为冠状动脉血管取样点，在冠状动脉血管取样点上添加冠状动脉血管的运动信息，以建立心动周期的七个阶段内的运动模型，得到它们在整个心脏周期的运动模型；然后对这些取样点在同一时刻的位置进行B样条曲线拟合得到对应于每一时刻的冠状动脉血管中轴三维静态模型，一个周期中所有时刻的冠状动脉血管中轴三维静态模型组成冠状动脉血管中轴三维动态模型；

(2.2)在冠状动脉血管中轴三维动态模型上加半径信息，得到冠状动脉血管三维动态模型；

第3步利用从不同个体的单臂X射线造影图像重建得到的三维冠状动脉血管树，补充和更新冠状动脉血管三维动态模型，其过程为：

(3.1)利用单臂X射线造影图像，构建造影***世界中心坐标系下的不同个体的三维冠状动脉血管树；在重建得到的三维冠状动脉血管树中找到与冠状动脉血管中轴三维静态模型中血管对应的血管段为参考血管段；

(3.2)对重建的三维冠状动脉血管树进行空间变换，计算变换后的参考血管段与其在冠状动脉血管中轴三维静态模型中的对应血管段的平均距离，当该距离最小时得到最佳变换矩阵；

(3.3)利用最佳变换矩阵将冠状动脉血管模型归一化到造影***世界中心坐标系后，计算重建的三维冠状动脉血管树的坐标矢量角，当坐标矢量角大于坐标矢量角门限h₁，或平均距离大于平均距离门限h₂时，则将该个体重建的三维冠状动脉血管树作为具有特异性的冠状动脉血管中轴三维静态模型，然后再加入半径信息和运动信息，得到具有特异性的冠状动脉血管三维动态模型；当坐标矢量角小于坐标矢量角门限h₁，且平均距离小于平均距离门限h₂，则将重建的三维冠状动脉血管树与冠状动脉血管中轴三维静态模型进行融合，得到更加通用化的冠状动脉血管中轴三维静态模型，然后再加入半径信息和运动信息，得到更加通用化的冠状动脉血管三维动态模型。

本发明方法通过人体的心脏切片图像，利用B样条曲面拟合方法建立心脏与冠状动脉血管的三维静态模型，然后通过心脏上每个点的法向和切向的运动建立冠状动脉血管三维动态模型，并利用从真实的单臂X射线造影图像重建得到的较详细冠状动脉血管结构补充或建立新的冠状动脉血管模型。具体而言，本发明具有如下四个方面的技术效果：

(1)通过心脏上各点的切向的旋转运动，获得心脏整体的旋转和局部扭转运动信息，进而得到附着于心脏表面的冠状动脉血管的运动信息，使得冠状动脉血管三维动态模型模型更接近现实情况；

(2)通过将建立的冠状动脉血管模型的运动与文献中的做标记实验测得结果进行比较，验证我们所建立冠状动脉血管模型的正确性，并在两者出现较大差异时，能以相对真实的文献数据作为参考来指导冠状动脉血管模型的递推建模，使冠状动脉血管模型的运动不断趋近于现实中人体冠状动脉血管的运动；

(3)在利用一组人体的心脏切片图像建立的冠状动脉血管模型基础上，根据实际单臂X射线造影图重建个体的三维冠状动脉血管树，并通过图像分析上述三维冠状动脉血管树与冠状动脉血管模型结构的差异，补充原有冠状动脉血管模型细节拓扑结构，并可通过拓扑结构差异分析建立具有特异性的冠状动脉血管模型，对于冠状动脉血管的共性与个体差异性都有很好的诠释；

(4)通过建立人体心脏不同年龄段的动态(或静态)的冠状动脉血管模型，分析人体心脏在不同年龄段发生的变化，逐步建立不同人种、不同个体、不同年龄段的具有特异性的人体动态心脏模型，对冠状动脉血管乃至心血管的发展和病变进行分析，为我们的临床治疗或研究提供客观的指导和参考，从而实现定量的医学诊断。

附图说明

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

图2是建立心脏三维动态模型流程图；

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

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

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

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

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

图5是心动周期各时期图

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

图7是文献中心脏心底和心尖部分的旋转角度随时间变化的曲线；

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

图11是左冠状动脉主支的根节点(LM0)，左旋支和前室间支的交叉点(LMb)，左旋支上离根部5mm的血管点(LCX5)，前室间支上离根部5mm的血管点(LAD5)四个点分别在方向后前(Posterior-Anterior)，下上(Inferior-Superior)，和左右(Left-Right)上一个心动周期内的运动曲线；

表1是图11中四个特征点的运动位移；

表2是文献中被选取四个对应于表1的特征点的运动位移；

图12是冠状动脉血管三维动态模型更新流程图；

图13(a)是个体A重建的三维血管树结构图；

图13(b)是个体B重建的三维血管树结构图；

图13(c)是个体C重建的三维血管树结构图；

图13(d)是个体D重建的三维血管树结构图；

图13(e)是冠状动脉血管模型三维结构图；

图14(a)是个体A与模型的三维融合图；

图14(b)是在图13(a)的基础上，进行细节补充，其中红色血管段即为添加血管。

具体实施方式

本发明方法通过同时考虑心脏的法向和切向的运动，在建立心脏模型的基础上构建冠状动脉血管模型，并利用不同个体的单臂X射线造影图像数据，重建个体冠状动脉血管的三维树，然后将重建后的结构与冠状动脉血管模型作差异比较，若拓扑差异较大，则构建新的具有特异性的冠状动脉血管模型，反之，则对两者进行融合，以获得更具一般性的冠状动脉血管模型。

心脏的运动是十分复杂的，根据医学观察和生物医学工程的实验结果证明，室壁运动由心肌的内向收缩运动、心脏的水平移动及心脏绕轴的转动三种主要运动方式组成，此外还有局部扭转、局部伸展等局部运动。由于心脏的水平移动不会造成心脏的形变，同时其运动也较小，这里暂对其不予考虑。在建立心脏三维动态模型时，我们主要考虑心脏的收缩与旋转两种运动。

冠状动脉血管三维骨架模型(即冠状动脉血管中轴三维静态模型)加入半径信息构成冠状动脉血管三维静态模型；冠状动脉血管中轴三维静态模型加入运动信息后构成冠状动脉血管中轴三维动态模型(即冠状动脉血管中轴的周期运动模型)，在冠状动脉血管中轴三维动态模型加入半径信息构成冠状动脉血管三维动态模型，上述各类关于冠状动脉血管的模型统称为冠状动脉血管模型。

在建模前我们先对心脏的运动做如下合理假设：

①由于心脏的心底部分与大血管相连，因此可将心底看作是固定的；

②左右心房和左右心室的运动分别同步；

③心房与心室之间的运动不会互相影响；

④心包随着腔室运动；

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

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

⑦心脏腔室和心包表面上的各点的运动由法向和切向两个方向的运动组成；

⑧设心动周期为0.8秒，分别按下述八个时间点0，0.1s，0.15s，0.25s，0.40s，0.47s，0.58s，0.8s(记为t_e，e∈[0，7])将心动周期分为7个时间段，并合理假定每一时间段内心脏腔室和心包都近似为线性运动。

其中，①、③和⑦是在CN101301207A的基础上，针对于心脏切向运动的添加而提出的。

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

本发明方法通过人体心脏切片建立心脏三维静态模型与冠状动脉血管模型，然后根据心脏上每个点的法向和切向的运动来建立心脏运动模型，并根据冠状动脉血管与心脏的附着关系，建立动态的冠状动脉血管模型。同时，本发明还根据从实际的X射线造影图中重建得到的三维冠状动脉血管树补充或建立新的具有个体特异性的冠状动脉血管模型。

以下结合附图，并以左冠状动脉血管为例对本发明方法的实现过程作进一步说明：

如图1所示，本发明方法包括以下几个步骤：

(1)建立心脏三维静态模型

通过对提取的每层心脏切片轮廓采样，获取心脏表面的三维采样点，再利用B样条对这些采样点进行曲面拟合，建立心脏三维静态模型，包括腔室静态模型和心包静态模型。

图3是将心脏腔室与表面合在一起构成心脏的三维静态模型。最外面的轮廓是心包，区域B1是左心室，区域B2是右心室，区域B4是左心房，区域B3是右心房。图4是通过改变控制点的位置模拟的具有个体特异性的心脏。

(2)构建腔室动态模型

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

腔室上的各点运动由法向和切向两个方向的运动组成，即 $\stackrel{\→}{D}\left(t\right)=\stackrel{\→}{r}\left(t\right)+\stackrel{\→}{m}\left(t\right),$ 其中表示腔室上的任意一点的整体运动矢量，表示该点的法向运动矢量，

表示该点的切向运动矢量，t表示时间。

下面分别讲述这两种运动矢量的求取方法：

(a)法向运动矢量

的求取

对于法向的运动，其运动方向就是腔室表面上的点的法线方向，当腔室收缩时，沿法线反方向向里运动，当其舒张时，沿法线方向向外运动。其运动大小是通过腔室体积的变化推导得到的腔室的形变，即 $V\left(t\right)\⇒r_v\left(t\right),$ 这里V(t)表示时间t时腔室的体积，r_v(t)表示腔室在法向方向的形变。具体步骤为：

1)在时间t时，根据曲线V(t)求出相对于时间t_pre的体积变化ΔV，t_pre表示原始的心脏三维静态模型在心动周期中对应的时间，t_pre在t₀与t₇之间取值。

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

3)计算从时间t_pre到t时每个点产生的沿着该点法线方向的位移量为d＝ΔV/S。

4)对曲面上(不包括腔室连接处)的每一个取样点P₁，计算在点P₁处的曲面法线方向

那么，

则为该点对应时间t下的法向运动矢量

在左心室中，时间t时每个点产生的沿着该点法线方向的位移量d可通过如下公式计算得出：

$d=\left\{\begin{array}{cc}210t/S\left(t\right),& t\∈[t_0,t_1)\\ 0,& t\∈[t_1,t_2)\\ 470(0.15-t)/S\left(t\right),& t\∈[t_2,t_3)\\ \frac{460}{3}(0.25-t)/S\left(t\right),& t\∈[t_3,t_4)\\ 0,& t\∈[t_4,t_5)\\ \frac{4600}{11}(t-0.47)/S\left(t\right),& t\∈[t_5,t_6)\\ \frac{150}{11}(t-0.58)/S\left(t\right)& t\∈[t_6,t_7)\end{array}\right.,$

式中，S(t)为时间t时左心室表面(不包括腔室连接处)的曲面面积。

(b)切向运动矢量

的求取

对于切向的运动，它是绕着心脏的中心轴(实际建模中，一般可看作为左心室轴)进行旋转运动的。设心脏的中心轴为l，它的确定可通过以下几步实现：

1)在已构建好的左心室静态模型上找到心尖处中心点A，并计算出左心室静态模型的体心B，连接A、B，做

的延长线，令其与左心室腔室壁的交点为C，则获得心脏的中心轴l的初始值矢量

2)在线段AC上等间隔取k个点，k为大于等于3的正整数，经过这k个点做k个垂直于的垂面，然后找到每个垂面中包含的左心室腔室壁点，记为p_i(p_i是k个垂面中的任意一个垂面与左心室腔室壁相交处的任意一个点，i是腔室壁点的序号)，再分别计算出点p_i与线段AC上对应取样点的距离，记为L_i，并求出L_i的方差σ_j(j是垂面的序号，j∈[1，k])，最后获得所有面的方差和 $E=\underset{j=1}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\σ}}_j;$

3)固定心尖点A，移动点B至B′，且要求

与

的夹角在0到

之间，再做

的延长线，令其与腔室壁的交点为C′，即获得矢量

那么同理，经2)中的运算，求出任意B′相对应的方差和E′；

4)当方差和取最小值，即E_min＝MIN(E，E′)时，E_min对应的中心轴即为心脏的中心轴l。

心脏的中心轴l确定后，心脏上某点P的旋转角度随时间t的变化为θ(t)(其旋转角度都是以舒张末期时点的位置作为参考)。则旋转运动之后点P的位置为 $\stackrel{\→}{{\mathrm{OP}}^{\′}}=\stackrel{\→}{\mathrm{OP}}*R_l\left(\mathrm{\θ}\right),$ 其中P′表示旋转后的点，R_l(θ)表示以l为旋转轴，旋转角度为θ的旋转矩阵，可通过下式进行计算：

R_l(θ)＝R_l→ZR_Z(θ)R_l→Z ^-1，

其中R_l→Z表示将直线l变到Z轴的旋转矩阵，R_Z(θ)表示绕Z轴的旋转θ角度的旋转矩阵，那么切向运动矢量 $\stackrel{\→}{m}\left(t\right)=\stackrel{\→}{{\mathrm{OP}}^{\′}}-\stackrel{\→}{\mathrm{OP}}.$

根据文献中对心脏切向运动的生理研究，从心底到心包，心脏上的各点的旋转运动大小都不一样。假设在心脏收缩运动过程中从心底到心尖各点的旋转角度是线性变化的，参考文献中的结果可以得到心尖和心底这两部分的冠状动脉血管各点在整个心脏运动周期中每个时间段的旋转角度，然后通过线性插值得到心脏上每一个点的旋转角度。文献中心脏心底和心尖部分的旋转角度随时间变化的曲线如图7所示。

根据公式 $\stackrel{\→}{D}\left(t\right)=\stackrel{\→}{r}\left(t\right)+\stackrel{\→}{m}\left(t\right)$ 得到腔室上的任意一点在心动周期八个时间点的整体运动矢量

(e为0到7间的整数)，利用各点的整体运动矢量

获得对应时间点t_e下的各点新位置，即得到t_e时的腔室静态模型，那么，将时间点t₀到时间点t₇的腔室静态模型采用线性关系连接起来，从而得到腔室动态模型。

(3)构建心包动态模型：

心包的建模是基于腔室动态模型和心肌收缩模型，其中，心肌收缩模型我们定义为，在收缩过程中心肌体积保持不变。令心脏在收缩前、后的体积为V₁和V₂，H₁和H₂是心肌收缩前、后在心包和腔室之间的厚度，L′₁和L′₂是心脏在收缩前、后的心肌纤维长度。根据之前的假设，我们可以得到下列公式：

L′₁H₁ ²=L′₂H₂ ² $\frac{H_1}{H_2}={\left(\frac{V_1}{V_2}\right)}^{1/6}$

$\frac{V_1}{V_2}=\frac{L_1^{\′3}}{L_2^{\′3}}$

构建心包动态模型的具体步骤如下：

(3.1)从时间t_pre开始。t_pre表示原始的心脏三维静态模型在心动周期中对应的时间，t_pre在t₀与t₇之间取值。

(3.2)对心包上的每一个取样点P₂，计算在这一点的法线

找到法线

与腔室的交叉点Q。

(3.3)在时间t＝t_e时，根据腔室动态建模中求取法向运动矢量和切向运动矢量的方法得到交叉点Q的法向运动矢量

(与交叉点Q的法向运动矢量

相同)以及该时刻下的旋转角度θ(t_e)，t_e表示心动周期八个时间点中的任何一个，e∈[0，7]。

(3.4)根据心包和腔室之间的心肌层厚度计算心脏壁的厚度H＝|P₂Q |，然后得到新的心脏壁厚度 $H^{\′}=H{\left(\frac{V_e}{V_{\mathrm{pre}}}\right)}^{1/6},$ 其中，V_e为时间点t_e时的心脏体积，V_pre为时间t_pre时的心脏体积。

(3.5)令其中P′₂为形变之后的取样点P₂的新位置，O为模型所处坐标系的原点，R_l(θ(t_e))表示以l为心脏旋转轴，旋转角度为θ(t_e)的旋转矩阵，θ(t_e)为时间点t_e时的旋转角度。

(3.6)根据新的取样点构建心包模型，保存该心包模型及其对应的时间点t_e。

(3.7)移动到下一个时间点t_e+1，重复步骤(3.2)到步骤(3.6)，直到一个心动周期结束。

(3.8)假定相邻时间点间的心包上的点的运动为线性运动，将对应于选取时间点t_e的心包模型连接起来，构成心包动态模型。

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

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

(4)建立冠状动脉血管三维静态模型

建立冠状动脉三维静态模型是为了建立冠状动脉的三维运动模型，其步骤是：

(4.1)提取冠状动脉血管切片轮廓；

根据解剖学的知识和人体横断解剖的图片，在每幅原始图片中分出冠状动脉血管。图9是原始心脏切片图像及其冠状动脉血管提取结果对比图。

(4.2)提取冠状动脉血管的中轴点与半径；

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

(4.3)重建冠状动脉血管三维骨架模型：把步骤(4.2)提取的冠状动脉血管各个分支的中轴点用B样条曲线进行拟合得到冠状动脉血管三维骨架模型，即初始时刻下的冠状动脉血管中轴三维静态模型。

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

(5)建立冠状动脉血管三维动态模型：

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

(5.1)建立冠状动脉血管中轴三维动态模型；

在冠状动脉血管中轴上面取等间隔的点作为冠状动脉血管取样点，由于冠状动脉血管附着于心包表面，因此，可按照心包运动的规律，即将心包上距离某上述取样点最近的点的运动作为该取样点的运动，以获得冠状动脉血管取样点在不同时间的新位置。根据心脏生理学知识，心动周期分为7个阶段，合理地假设这些取样点在这7个阶段里分别近似为线性运动。对冠状动脉血管取样点通过分别建立这7个阶段内的动态模型，得到它们在整个心动周期的的动态模型。然后对这些取样点在同一时刻的位置进行B样条曲线拟合得到对应于每一时刻的冠状动脉血管中轴三维静态模型，一个周期中所有时刻的冠状动脉血管中轴三维静态模型就组成了冠状动脉血管中轴三维动态模型。图10是建立的一组对应七个阶段的冠状动脉血管中轴三维动态模型。

(5.2)在冠状动脉血管中轴三维动态模型上加半径信息

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

为了对冠状动脉血管三维动态模型的运动作进一步分析，可以对冠状动脉血管三维动态模型进行验证，其方法为：

在左冠状动脉血管模型上选取4个点：左冠状动脉主支的根节点(LM0)，左旋支和前室间支的交叉点(LMb)，左旋支上离根部5mm的血管点(LCX5)，前室间支上离根部5mm的血管点(LAD5)，在生理坐标系下的后前(Posterior-Anterior)，下上(Inferior-Superior)，和左右(Left-Right)三个方向分析其在一个心动周期中的运动，并与文献中相应四点在这三个方向的位移数据进行比较。

图11所示分别为动态的左冠状动脉血管模型验证中这四个点在一个周期内的运动曲线，表1中所示为模型中这些点在三个方向上以及三维空间中的最大位移，表2中所示为文献中的结果，该结果是通过对10个病人的双臂X射线造影图序列进行标记后提取获得的，具有较高的可信度。

由两表中数据的对比可知，该左冠状动脉血管模型上血管点的运动规律与文献中标记点的运动是一致的。除了在下上方向的位移较文献结果略大些外，其余数据值均落在它们相应的置信区间内。因此，总体上来说，该左冠状动脉运动模型是可接受的，它符合左冠状动脉运动的一般规律。

(6)利用从不同个体的单臂X射线造影图像重建得到的三维冠状动脉血管树，补充和更新冠状动脉血管三维动态模型。

冠状动脉血管既存在共性也存在个体差异性，我们可以根据不同个体的共性按照上面的方法建立动态的冠状动脉血管模型，但是由于每个个体还是存在一定的差异性，也需要建立一些具有个体特异性的冠状动脉血管模型，因此可以通过对从实际的单臂X射线造影图像重建得到的三维冠状动脉血管树的分析，当三维冠状动脉血管树表现出的个体特异性较大时，可作为具有个体特异性的冠状动脉血管模型。

在重建出三维冠状动脉血管树后，将三维冠状动脉血管树与冠状动脉血管模型进行差异比较，若差异较大，则构建新的具有特异性的冠状动脉血管模型，反之，则对两者进行融合，以获得更具一般性的冠状动脉血管模型，并补充冠状动脉血管模型中不存在的血管。具体按照如下步骤进行：

(6.1)利用单臂X射线造影图像，构建造影***世界中心坐标系下的不同个体的三维冠状动脉血管树；在重建得到的三维冠状动脉血管树中找到与冠状动脉血管中轴三维静态模型中血管对应的血管段为参考血管段。由于在用模型指导造影图中冠状动脉血管标记时得到过造影图中冠状动脉血管与模型中冠状动脉血管的对应关系，三维冠状动脉血管中轴三维静态模型上的血管段与重建得到的三维冠状动脉血管树上血管段的对应关系也因此可以确定；

(6.2)对重建的三维冠状动脉血管树进行空间变换(包括旋转、平移)，计算变换后的参考血管段与其在冠状动脉血管中轴三维静态模型中的对应血管段的平均距离，当该距离最小时得到最佳变换矩阵；

(6.3)将重建的三维冠状动脉血管树的坐标矢量角(该参数与左冠状动脉树中左旋支和前室间支的交叉点LMb或者右冠状动脉树中锐缘支和后降支的交叉点RMb相关，反映了在生理坐标系中左冠状动脉树或者右冠状动脉树相对于心脏的分布情况，)以及上步计算得到的对应血管段的平均距离作为差异性判据。

上述坐标矢量角门限和平均距离门限的获取过程如下：

(6.3.1)利用最佳变换矩阵，将冠状动脉血管中轴三维静态模型归一化到造影***世界中心坐标系下，分别计算冠状动脉血管中轴三维静态模型和重建的三维冠状动脉血管树在此坐标系下的LMb或RMb，记为M，I_g，其中，点M为冠状动脉血管中轴三维静态模型LMb或RMb，点I_g为不同个体重建的三维冠状动脉血管树的LMb或RMb，g为不同个体对应的编号；

(6.3.2)记造影***世界中心坐标系的原点为O₁，计算矢量

与矢量

的夹角

即坐标矢量角，然后求和取平均，即

结果则作为坐标矢量角门限h₁，n是不同个体的个数；

(6.3.3)计算不同个体重建的三维冠状动脉血管树与冠状动脉血管中轴三维静态模型间对应血管段的平均距离，然后求和再取平均，那么，再取平均所得到的结果即为平均距离门限h₂；

在实际操作中，坐标矢量角门限h₁可以取值为11.1°，平均距离h₂可以取值为14.6。

(6.4)差异较大，也就是当求得的坐标矢量角大于坐标矢量角门限h₁，或距离大于平均距离门限h₂时，则将该个体重建的三维冠状动脉血管树作为具有特异性的冠状动脉血管中轴三维静态模型；差异不大时，即满足求得的坐标矢量角小于坐标矢量角门限h₁，且距离小于平均距离门限h₂，则将重建的三维冠状动脉血管树与冠状动脉血管中轴三维静态模型进行融合，得到更加通用化的冠状动脉血管中轴三维静态模型，其过程如下：

(6.4.1)将通过最佳变换矩阵变换后的参考血管段与其在冠状动脉血管中轴三维静态模型中的对应血管段以相同数目进行等间隔取样，并使这些取样点一一对应，将对应点对的坐标取平均得到一系列新的取样点，对这些点进行B样条曲线拟合得到新的冠状动脉血管中轴三维静态模型；

(6.4.2)对冠状动脉血管中轴三维静态模型细节进行补充。当出现冠状动脉血管中轴三维静态模型中不存在的血管结构时，如果它在对不同个体的单臂X射线造影图进行重建的三维冠状动脉血管树中都存在，则将它添加到冠状动脉血管中轴三维静态模型中。在添加的过程中，保持血管的相对走向不变，将整段血管连接到融合后的新冠状动脉血管中轴三维静态模型上，从而对冠状动脉血管中轴三维静态模型进行补充。

(6.5)在步骤(6.4)得到的具有特异性或通用化的冠状动脉血管中轴三维静态模型中加入半径信息和运动信息，得到相应的具有特异性或通用化的冠状动脉血管三维动态模型。

图12为冠状动脉血管三维动态模型更新的流程图；

利用最佳变换矩阵将冠状动脉血管模型归一化到造影***世界中心坐标系后，计算重建的三维冠状动脉血管树的坐标矢量角，当坐标矢量角大于坐标矢量角门限h₁，或平均距离大于平均距离门限h₂时，则将该个体重建的三维冠状动脉血管树作为具有特异性的冠状动脉血管中轴三维静态模型，然后再加入半径信息和运动信息，得到冠状动脉血管三维动态模型；当坐标矢量角小于坐标矢量角门限h₁，且平均距离小于平均距离门限h₂，则将重建的三维冠状动脉血管树与冠状动脉血管中轴三维静态模型进行融合，得到更加通用化的冠状动脉血管中轴三维静态模型，然后再加入半径信息和运动信息，得到冠状动脉血管三维动态模型。

通过四个不同个体(分别记为A、B、C和D)重建的三维冠状动脉血管树与冠状动脉血管模型的分析，发现冠状动脉血管模型与个体A和C重建的三维冠状动脉血管树基本的较大尺度的较高级别结构很相似，只是重建得到的血管结构分支即较小尺度较低级别的结构比冠状动脉血管模型中要多，因此可以通过重建的三维冠状动脉血管树对冠状动脉血管模型进行补充。另外，个体B重建得到的三维冠状动脉血管树与冠状动脉血管模型相差较大，表现在左旋支较短，且分支较少，而由个体D重建的三维冠状动脉血管树则显示该个体的心包很狭窄，因而可以分别将C、D分别作为新的具有个体特异性的冠状动脉血管模型，加入到模型库中。

图13中分别为个体A，B，C，D重建的三维冠状动脉血管树及冠状动脉血管模型三维结构图。

图14中分别为个体A与冠状动脉血管模型的融合图及其细节补充图(虚线血管即为添加血管)。

以上所述为本发明的较佳实施例而已，但本发明不应该局限于该实施例和附图所公开的内容。所以凡是在不脱离本发明所公开的精神下完成的等效或修改，都落入本发明保护的范围。

Claims (7)

1.一种人体冠状动脉血管的递推建模方法，其特征在于，该方法包括下述步骤：

第1步建立冠状动脉血管三维静态模型，其过程为：

(1.1)提取冠状动脉血管切片轮廓；

(1.2)提取冠状动脉血管的中轴点与半径，中轴点为冠状动脉血管切片轮廓的形心，同时假设冠状动脉为圆形，冠状动脉血管的半径为轮廓的内切圆半径；

(1.3)重建冠状动脉血管三维骨架模型：把步骤(1.2)提取的冠状动脉血管各个分支的中轴点用B样条曲线进行拟合得到冠状动脉血管三维骨架模型，即初始时刻下的冠状动脉血管中轴三维静态模型；

(1.4)构建带半径的冠状动脉三维模型：在所述的冠状动脉血管三维骨架模型上使用广义圆柱模型加入半径信息，建立带半径的冠状动脉三维模型，即冠状动脉血管三维静态模型；

第2步建立冠状动脉血管三维动态模型，其过程为：

(2.1)建立冠状动脉血管中轴三维动态模型：

利用冠状动脉血管与心包的附着关系，将心包动态模型的运动信息转化为冠状动脉血管的运动信息，该运动信息包括法向运动矢量和切向运动矢量；在冠状动脉血管中轴取等间隔的点作为冠状动脉血管取样点，在冠状动脉血管取样点上添加冠状动脉血管的运动信息，以建立心动周期的七个阶段内的运动模型，得到它们在整个心脏周期的运动模型；然后对这些取样点在同一时刻的位置进行B样条曲线拟合得到对应于每一时刻的冠状动脉血管中轴三维静态模型，一个周期中所有时刻的冠状动脉血管中轴三维静态模型组成冠状动脉血管中轴三维动态模型；

(2.2)在冠状动脉血管中轴三维动态模型上加半径信息，得到冠状动脉血管三维动态模型；

第3步利用从不同个体的单臂X射线造影图像重建得到的三维冠状动脉血管树，补充和更新冠状动脉血管三维动态模型，其过程为：

(3.1)利用单臂X射线造影图像，构建造影***世界中心坐标系下的不同个体的三维冠状动脉血管树；在重建得到的三维冠状动脉血管树中找到与冠状动脉血管中轴三维静态模型中血管对应的血管段为参考血管段； 

(3.2)对重建的三维冠状动脉血管树进行空间变换，计算变换后的参考血管段与其在冠状动脉血管中轴三维静态模型中的对应血管段的平均距离，当该距离最小时得到最佳变换矩阵；

(3.3)利用最佳变换矩阵将冠状动脉血管模型归一化到造影***世界中心坐标系后，计算重建的三维冠状动脉血管树的坐标矢量角，当坐标矢量角大于坐标矢量角门限h₁，或平均距离大于平均距离门限h₂时，则将该个体重建的三维冠状动脉血管树作为具有特异性的冠状动脉血管中轴三维静态模型，然后再加入半径信息和运动信息，得到具有特异性的冠状动脉血管三维动态模型；当坐标矢量角小于坐标矢量角门限h₁，且平均距离小于平均距离门限h₂，则将重建的三维冠状动脉血管树与冠状动脉血管中轴三维静态模型进行融合，得到更加通用化的冠状动脉血管中轴三维静态模型，然后再加入半径信息和运动信息，得到更加通用化的冠状动脉血管三维动态模型。

2.根据权利要求1所述的人体冠状动脉血管的递推建模方法，其特征在于，步骤(2.1)中按照下述过程建立心包动态模型：

(2.1.1)从时间t_pre开始；t_pre表示原始的心脏三维静态模型在心动周期中对应的时间；

(2.1.2)对心包的每一个取样点P₂，计算在这一点的法线 

找到法线 

与腔室的交叉点Q；

(2.1.3)在时间t＝t_e时，在腔室动态模型中求取交叉点Q在时间点t_e的法向运动矢量 

以及时间点t_e时的旋转角度θ(t_e)，t_e表示心动周期八个时间点中的任何一个，e∈[0，7]；

(2.1.4)根据心包和腔室之间的心肌层厚度计算心脏壁的厚度H＝|P₂Q|，然后得到新的心脏壁厚度 

其中，V_e为时间点t_e时的心脏体积，V_pre为时间t_pre时的心脏体积；

(2.1.5)令其中P′₂为形变之后的取样点P₂的新位置，O为模型所处坐标系的原点，R_l(θ(t_e))表示以l为心脏旋转轴，旋转角度为θ(t_e)的旋转矩阵，θ(t_e)为时间点t_e时的旋转角度； 

(2.1.6)以取样点P₂的新位置作为新的取样点，构建心包模型，保存该心包模型及其对应的时间点t_e；

(2.1.7)移动到下一个时间点t_e+1，重复步骤(2.1.2)到步骤(2.1.6)，直到一个心动周期结束；

(2.1.8)假定相邻时间点间的心包上的点的运动为线性运动，将对应于选取时间点t_e的心包模型连接起来，构成心包动态模型。

3.根据权利要求2所述的人体冠状动脉血管的递推建模方法，其特征在于，步骤(2.1.3)中腔室动态模型的建立过程为：

(A1)求取腔室上的各点的法向运动矢量 

和切向运动矢量 

(A2)根据公式 

得到腔室上的任意一点的整体运动矢量 

利用各点的整体运动矢量 

得到腔室动态模型。

4.根据权利要求3所述的人体冠状动脉血管的递推建模方法，其特征在于，步骤(2.1)按照下述过程求取法向运动矢量 

(B1)在时间t时，计算腔室在时间t相对于时间t_pre的体积变化ΔV，t_pre表示原始心脏三维静态模型在心动周期中对应的时间；

(B2)计算不包括腔室连接处在内的左心室表面的曲面面积S；

(B3)计算从时间t_pre到t后每个点产生的沿着该点法线方向的位移量为d＝ΔV/S；

(B4)对不包括腔室连接处在内的左心室表面上的每一个取样点P₁，计算在点P₁处的曲面法线方向 

则为该点对应时间t下的法向运动矢量 

5.根据权利要求3所述的人体冠状动脉血管的递推建模方法，其特征在于，步骤(2.1)按照下述过程求取切向运动矢量 

(C1)在左心室静态模型上找到心尖处中心点A，并计算出左心室静态模型的体心B，连接A、B，做 

的延长线，令其与左心室腔室壁的交点为C，则获得心脏的中心 轴l的初始值矢量 

(C2)在线段AC上等间隔取k个点，k为大于等于3的正整数，经过这k个点做k个垂直于 的垂面，然后找到每个垂面中包含的左心室腔室壁点，记为p_i，i是腔室壁点的序号，再分别计算出点p_i与线段AC上对应取样点的距离，记为L_i，并求出L_i的方差σ_j，j是垂面的序号，j∈[1，k]，获得所有垂面的方差和 

(C3)固定心尖点A，移动点B至B′，且要求 与 

的夹角在0到 之间，再做 

的延长线，令其与腔室壁的交点为C′，即获得矢量 再利用步骤(C2)的方法计算出任意B′相对应的方差和E′；

(C4)当方差和取最小值，即E_min＝MIN(E，E′)时，E_min对应的中心轴即为心脏的中心轴l；

(C5)心脏的中心轴l确定后，心脏上任一点P的旋转角度随时间t的变化为θ(t)，设旋转运动后的点置为P′，模型所处坐标系的原点为O，点P′的位置为 R_l(θ)表示以l为旋转轴、旋转角度为θ的旋转矩阵，切向运动矢量 

6.根据权利要求1至5中任一所述的人体冠状动脉血管的递推建模方法，其特征在于，第3.3步中，坐标矢量角门限h₁和平均距离门限h₂的获取过程为：

(D1)将冠状动脉血管中轴三维静态模型归一化到造影***世界中心坐标系下，分别计算冠状动脉血管中轴三维静态模型和重建的三维冠状动脉血管树在此坐标系下两主支的交叉点，记为M，I_g，其中，点M为左冠状动脉血管模型左旋支和前室间支的交叉点LMb，或者为右冠状动脉血管模型锐缘支和后降支的交叉点RMb，点I_g为不同个体重建的三维冠状动脉血管树的LMb或RMb，g为不同个体对应的编号；

(D2)记造影***世界中心坐标系的原点为O₁，计算矢量 

与矢量 的夹角 

即坐标矢量角，然后求和取平均，即 

结果则作为坐标 矢量角门限h₁，n是不同个体的个数；

(D3)计算不同个体重建的三维冠状动脉血管树与冠状动脉血管中轴三维静态模型间对应血管段的平均距离，然后求和再取平均，那么，再取平均所得到的结果即为平均距离门限h₂。

7.根据权利要求1至5中任一所述的人体冠状动脉血管的递推建模方法，其特征在于，第3.3步中，重建的三维冠状动脉血管树与冠状动脉血管中轴三维静态模型进行融合的过程如下：

(E1)将通过最佳变换矩阵变换后的三维冠状动脉血管树的参考血管段与其在冠状动脉血管中轴三维静态模型中的对应血管段以相同数目进行等间隔取样，并使这些取样点一一对应，将对应点对的坐标取平均得到一系列新的取样点，对这些点进行B样条曲线拟合得到新的冠状动脉血管中轴三维静态模型；

(E2)对冠状动脉血管中轴三维静态模型细节进行补充；当出现冠状动脉血管中轴三维静态模型中不存在的血管结构时，如果它在对不同个体的单臂X射线造影图进行重建的三维冠状动脉血管树中都存在，则将它添加到冠状动脉血管中轴三维静态模型中；在添加的过程中，保持血管的相对走向不变，将整段血管连接到融合后的新冠状动脉血管模型上，从而对冠状动脉血管中轴三维静态模型进行补充。 

Images