专利名称：一种冠状动脉血管轴线的四维重建方法
技术领域：
本发明涉及一种根据覆盖一个或多个心动周期的数字X射线冠状动脉造影图像序列对血管轴线进行四维(三维+时间)重建的方法，属医学检测技术领域。
背景技术：
X射线冠状动脉造影(X-ray coronary angiography，CAG)是目前临床广泛采用的诊断和治疗冠心病的影像手段。CAG诊断的最大特点在于可视性，即通过静止或动态观察造影剂的充盈和消失情况来判断心血管解剖学形态异常的部位、性质和程度。X射线造影成像是将空间血管结构重叠到二维成像平面上，即投影成像，从而丢失了临床诊断中所需的大部分三维空间信息。传统的诊断过程中，医生需要利用解剖、病理等专门知识和临床经验，根据血管在多个方向的投影想象其三维形态，从而做出评价。分析结果的准确程度密切依赖于医生的临床经验和专业知识，客观性和可重复性差。同时，冠状动脉附着在心外膜表面上，随着心脏有节律地收缩和舒张，因而在心动周期的不同时刻，血管的形态结构和空间位置会发生很大变化。单纯依靠一个时刻的二维投影或三维重建结果很难得到对血管形态结构的准确描述。
在一个或多个心动周期中拍摄的X射线冠状动脉造影图像序列记录了心脏搏动时的血管变形和心脏供血功能等动态信息。采用数字图像处理的方法，可以从不同时刻、不同角度拍摄的图像序列中客观、定量地提取隐含的空间和时间信息，并以恰当的方式表达出来。这些信息对于冠心病的临床诊治具有重要的辅助作用。同时冠脉的四维(三维+时间)重建也是造影图像与其它冠脉成像技术(如血管内超声(IVUS)、血管内OCT等)的融合、心动周期中的冠脉运动估计等的关键步骤。
为了实现对冠脉在整个心动周期中的四维重建，现有的方法多是按照“自底向上”，即2D→3D的顺序，逐个时刻进行三维重建。首先从两个角度的造影图像中提取出主要血管分支的轴线；之后，进行左右两个角度间轴线对应点的匹配；最后，根据两个投影点的坐标，计算出各轴线点的三维坐标。由于需要首先从原始图像中提取出单像素、八连通的血管轴线，因此结果的精度在很大程度上取决于二维提取的准确度。但临床采集的造影图像通常有较严重的噪声污染，而且图像中经常会出现血管重叠的现象，这种情况下仅从一个角度的图像很难得到准确的血管轮廓。目前全自动的提取算法尚不能对何种质量的图像均保成功，仍需操作者的手动参与。同时，两个角度之间血管轴线投影点的匹配直接影响到重建的精度。对于离散的二维轴线，一般是利用立体成像中的外极线约束逐点匹配，该方法很难达到较高的精度，重建结果的连续性不好，计算量也很大。此外，由于需要对序列中每个时刻的图像均重复上述步骤，因此要完成整个序列的血管跟踪，所需的计算成本过高，且需要操作者大量的手动参与，实用价值较低。
发明内容
本发明的目的是克服现有技术的不足、提供一种操作简便、计算量小、工作效率高的冠状动脉血管轴线的四维重建方法。
本发明所称问题是以下述技术方案实现的 一种冠状动脉血管轴线的四维重建方法，它是在采集两个角度的覆盖一个或多个心动周期的X射线冠状动脉造影图像序列后，首先建立X射线血管造影系统在两个角度的投影模型，并推导出两幅不同角度的造影图像之间的几何变换关系，再对手动选取的第一时刻图像中感兴趣血管段的采样点进行三维重建，并以连接各重建点所得折线作为snake模型的初始位置，通过snake变形获得该段血管第一时刻的三维轴线，而对于图像序列中的每一后续时刻，则均以其前一时刻的3D血管轴线作为snake模型的初始位置，通过snake变形获得三维血管轴线，从而完成整个序列血管轴线的四维重建，具体步骤如下 a、以心电信号为同步信号，采集两个角度的、覆盖一个或多个心动周期的单面X射线冠状动脉造影同步图像序列，并对原始图像进行初步的滤波、去噪、畸变校正等处理，增强视觉效果； b、建立X射线血管造影系统在两个角度的投影模型，推导两幅不同角度的造影图像之间的几何变换关系 设点s1和s2为两次造影过程中X射线源焦点的位置，分别以s1和s2为原点，建立空间坐标系s1x1y1z1和s2x2y2z2；坐标系U1V1O1和U2V2O2为成像平面坐标系；D1和D2分别是s1和s2到各自成像平面的垂直距离，空间血管上的点P在成像平面A和B上的投影点分别为p1(u1，v1)和p2(u2，v2)， 则s1x1y1z1和s2x2y2z2之间的几何变换关系为 其中，R3×3为旋转矩阵R＝RY(α2)·RX(β2)·RX(-β1)·RY(α1)
为平移向量 L1和L2为X射线源焦点s1和s2到旋转中心的距离；(α1，β1)和(α2，β2)分别是A和B的造影角度；(x1，y1，z1)和(x2，y2，z2)分别表示空间点P在坐标系s1x1y1z1和s2x2y2z2中的坐标； 根据透视投影成像的几何关系可知，空间点的三维坐标与其二维投影坐标之间的关系为 其中， c、在序列中第一时刻的图像中，通过人工取点获得感兴趣血管段在左右投影中的近似中心线，用折线表示，并利用外极约束得到两个角度间对应点的匹配； d、对手动选取的采样点进行三维重建； e、连接各3D点，所得折线作为3D snake的初始位置，通过使预先设定的能量函数最小，snake在空间中变形，最终得到具有最小能量的最优位置，就是第一时刻的3D血管轴线； f、后续时刻血管轴线的四维重建 对于序列中的后续时刻，以前一时刻的3D血管轴线作为当前时刻snake的初始位置，通过snake变形，完成整个序列中血管轴线的四维重建。
上述冠状动脉血管轴线的四维重建方法，所述X射线冠状动脉造影图像序列的两个采集角度之间夹角的取值范围为60°～120°。
上述冠状动脉血管轴线的四维重建方法，在第一时刻图像的感兴趣血管段中选取的三维重建采样点包括起点、终点和3～6个中间点。
本发明按照“自顶向下”的方式，采用3D snake模型技术，通过使能量函数最小，表示血管轴线的曲线直接在三维空间中变形，完成血管轴线的四维重建。该方法将序列图像中血管的二维中心线提取、三维重建和相邻时刻间的三维运动跟踪集成到一个框架中完成，提高了运算精度和速度。与“自底向上”的传统重建算法相比较，由于曲线直接在三维空间中变形，因此避免了逐点匹配，提高了重建精度。同时操作者的参与被减小到只需在序列的首帧中选择感兴趣血管段上的若干点。本发明不仅大大减少了由操作者参与所引入的不确定性和误差，提高了结果的可重复性，而且操作简便，工作效率高。
图1是本发明的实施步骤流程图； 图2是本发明的单面造影系统在两个角度的成像示意图； 图3是本发明的血管造影系统的造影角度示意图； 图4是采用本发明的左冠前降支(LAD)首时刻轴线重建的实施例图像；其中，图4(a)是首时刻的造影图像对和初始的3D snake；图4(b)是首时刻的造影图像对、重建出的血管轴线及其在两个成像平面上的二维投影； 图5是采用本发明的左冠前降支(LAD)后续时刻轴线序列重建的实施例图像，其中，图5(a)～(d)分别是序列中的第4、6、8、10个时刻的血管轴线重建结果及其在两个成像平面上的二维投影。
图中各符号为A、B、成像平面；s1、s2、两次造影过程中X射线源焦点的位置；s1x1y1z1、以s1为原点的空间坐标系；s2x2y2z2、以s2为原点的空间坐标系；OXYZ、以旋转中心为原点的空间坐标系；U1V1O1、成像平面A上的直角坐标系；U2V2O2、成像平面B上的直角坐标系；D1、s1到成像平面A的垂直距离；D2、s2到成像平面B的垂直距离；P、空间血管上的点；p1、P点在成像平面A上的投影；p2、P点在成像平面B上的投影；u1、p1在坐标系U1V1O1内的横坐标；v1、p1在坐标系U1V1O1内的纵坐标；u2、p2在坐标系U2V2O2内的横坐标；v1、p2在坐标系U2V2O2内的纵坐标；K1、K2、外极线；(α1，β1)、成像平面A的造影角度；(α2，β2)、成像平面B的造影角度；L1、X射线源s1到旋转中心的距离；L2、X射线源s2到旋转中心的距离。
具体实施例方式 下面结合附图和实例对本发明作进一步详细说明。
如附图1所示，本发明方法的步骤包括 (1)图像采集和预处理 采用C型臂单面X射线血管造影系统获得两个角度的、至少覆盖一个心动周期的冠状动脉造影图像序列，要求两个角度之间的夹角在60°至120°之间。记录成像系统参数(造影角度，X射线源到成像平面的距离)。采用同步记录的心电信号选取心动周期中相同相位的两个角度的单面图像对。
对原始图像进行必要的预处理，主要包括畸变校正、均衡对比度、去除噪声和图像增强等，提高图像的视觉效果，为后续处理奠定基础。
(2)建立单面造影系统在两个近似正交角度的成像模型 如附图2所示，点s1和s2表示两次造影过程中X射线源焦点的位置。分别以s1和s2为原点，建立空间坐标系s1x1y1z1和s2x2y2z2；坐标系U1V1O1和U2V2O2为成像平面坐标系；D1和D2分别是s1和s2到各自成像平面的垂直距离，随成像面的移动而改变；空间血管上的点P在成像平面A和B上的投影点分别为p1(u1，v1)和p2(u2，v2)。
从坐标系s1x1y1z1到s2x2y2z2的变换运动是三维空间中的刚体运动，根据刚体运动理论，推导出坐标系s1x1y1z1和s2x2y2z2之间的几何变换关系 其中R3×3为旋转矩阵 R＝RY(α2)·RX(β2)·RX(-β1)·RY(α1) (2)
为平移向量 式中，L1和L2为X射线源s1和s2到旋转中心的距离；(α1，β1)和(α2，β2)分别是图像A和B的造影角度(如附图3所示)。距离和角度值均可从造影机上得到。将称为几何变换矩阵。
根据透视投影成像的几何关系可知 其中 (x1，y1，z1)和(x2，y2，z2)分别表示空间点P在坐标系s1x1y1z1和s2x2y2z2中的坐标。因此根据式(1)、(4)和(5)，由点p1和p2的坐标可反算出点P的三维坐标。
(3)第一时刻血管轴线的三维重建 在跟踪过程开始前，需要确定snake模型的初始位置。本发明采用手动取点的方法，首先从一个角度的图像中手动选取感兴趣血管段上的若干采样点(一般选取血管段的起点、终点和3～6个中间点即可，具体数目根据血管段的长度决定)。然后采用外极约束得到各点在另一个角度图像上的对应点。如附图2所示，空间点P与X射线源焦点s1和s2构成外极平面，该平面与图像平面A和B相交形成左右两条外极线K1和K2。根据外极约束原理，点P在图像A上的投影p1在图像B中的对应点p2一定位于K2上；P在图像B上的投影p2在图像A中的对应点p1一定位于K1上。由于几何变换矩阵和其它计算的误差，匹配点可能不会准确地出现在对应的外极线上，而在外极线的邻域内。本发明方法在该邻域内搜索和外极线距离最近的点作为匹配点。由这几组对应点分别求出它们的三维坐标。用直线段连接这些3D点，所得折线作为3D snake的初始位置。
snake是一条可变形的参数曲线c(s)＝(x(s)，y(s)，z(s))，s∈
。将该曲线离散化为N个点Ci(Xi，yi，Zi)(i＝1，2，…，N)，模型能量函数的离散表达式为 其中，Eint是内部能量 Eint(i)＝α(d-|ci-ci-1|)+β|ci-1-2ci+ci+1|(7) 保证曲线变形过程中的连续性和光滑性。式中第一项约束使点平均分布，既能满足连续性的要求，又不会产生曲线收缩的现象。在每次迭代结束时，点间的平均距离d的值被更新；第二项是曲率，使曲线保持平滑。α和β是权值。
外部能量Eext是保证模型收敛的外部力，吸引表示血管轴线的空间曲线移动到这样的位置3D曲线在左右两个图像平面上的投影位于图像上灰度最小、且灰度梯度最小的区域，即血管的2D中心线。根据式(1)和(5)，空间点P在图像A和B上的投影坐标(u1，v1)和(u2，v2)都可用该点的三维坐标c＝(x1，y1，z1)来表示 [u1v1]T＝TL(c)，[u2v2]T＝TR(c) (8) 其中TL和TR都是以c为自变量的函数。Eext包含两部分，分别对应于两个角度的二维投影 式中IL(TL(ci))＝IL(u1i，v1i)和IR(TR(ci))＝IR(u2i，v2i)分别是左右图像点的灰度值； 和分别是左右图像点的灰度梯度。由于造影图像中，血管的灰度值比背景小，所以这里权重参数γ取正值。而且感兴趣的是血管中心线，不是边缘，所以λ也取正值。
(4)后续时刻血管轴线的四维重建 得到第一时刻血管段的三维轴线之后，对于其后的序列图像，根据血管形状不突变的特点，把前一时刻的三维血管轴线作为当前时刻snake初始位置，实现对整个序列中血管轴线的四维重建。
在这一步骤中，本发明方法在snake能量函数中嵌入相似度匹配，通过对相邻时刻之间的目标进行配准，实现对血管的跟踪。假设在足够短的时间间隔内，运动前后图像中血管的灰度不发生变化，那么得到第一时刻的血管轴线之后，在开始对后续图像进行运动跟踪时，外部能量函数采用如下的形式 式中ILt和IRt分别表示时刻t左右图像点的灰度值；ILt-1和IRt-1分别表示时刻t-1左右图像点的灰度值；Δw是邻域窗口的大小。(10)式中的前两项与(9)式相同，表示吸引曲线向着图像中灰度最小、灰度梯度最小的区域移动的图像力。第三项是使得运动前后(时刻t-1和t)血管上对应点的Δw×Δw邻域内的灰度变化最小的外部约束力。
对于求解snake能量函数最小化的问题，在计算机视觉发展的最初阶段，一般都是采用变分法，迭代求解E-L偏微分方程。该方法不仅计算复杂，运算时间长，而且很容易收敛于局部最小值。采用动态规划(“Using dynamic programming for solving variational problems invision，”IEEE Transactions on pattern analysis and machine intelligence.vol.12，no.9，pp.855-867，1990)可避免上述问题，由于迭代向着总能量减少的方向进行，当总能量不再变化时，则迭代终止，因此保证了结果的收敛性。但是该方法的缺点是计算速度慢，需要较大的存贮空间。针对上述问题，本发明采用Williams(“A fast algorithm for active contours andcurvature estimation，”Computer Vision，Graphics and Image Processing，vol.55，no.1，pp.14-26，1992)提出的贪婪算法完成对最优snake的搜索，该算法具有动态规划方法的所有优点，并且计算速度快，所需存贮空间小。
对于重建出的各时刻的三维血管轴线，根据式(8)即可算出其在左右两个成像平面上的二维投影，即二维血管轴线。因此本发明方法将二维提取、三维重建和运动跟踪集成到一个框架中完成，提高了运算精度和速度。
附图4和5是对PHILIPS Integris CV单面全数字血管造影机临床拍摄到的左冠造影图像序列的实验结果。采集速率为15帧/秒、图像大小为512×512(像素)、灰阶为256、像素尺寸为0.3mm。图像序列的拍摄角度分别为RAO30°CAUD24°和LAO46°CRAN21°。图4是第一时刻两个角度的图像对，在RAO图像的左前降支(LAD)上选择八个点，根据外极约束可找到其在LAO图像中的匹配点。计算出它们的三维坐标，用直线段连接这些3D点，得到用折线表示的3D snake的初始形状，如图4(a)所示。通过使能量函数((6)，(7)和(9)式)最小，snake不断变形，最后得到第一时刻的3D血管轴线。采用三次B样条曲线对这些离散3D点进行拟合，得到一条连续曲线表示的血管轴线，它在两个成像平面上的投影就是左右图像中的二维血管轴线，见图4(b)(在原始图像中用白色曲线表示)。
对于其后的序列图像，把t时刻中模型的停留位置作为t+1时刻模型的初始位置，通过使能量函数((6)，(7)和(10)式)最小，模型在内外力作用下变形，实现血管轴线序列的三维重建，即四维重建。图5是选取的第4、6、8、10时刻的跟踪结果，同样在原始图像上用白色曲线表示二维血管轴线。由此结果可见，当心动周期结束时，即第十时刻的血管形状与第一时刻的形状(图4(b))很相似，这也验证了“在心动周期结束时，冠状动脉血管会恢复到其在周期开始时的状态”的结论。
权利要求
1、一种冠状动脉血管轴线的四维重建方法，其特征是，首先采集两个角度的覆盖一个或多个心动周期的单面X射线冠状动脉造影图像序列、建立X射线血管造影系统在两个角度的投影模型，并推导出两幅不同角度的造影图像之间的几何变换关系，再对手动选取的第一时刻图像中感兴趣血管段的采样点进行三维重建，并以连接各重建点所得折线作为snake模型的初始位置，通过snake变形获得该段血管第一时刻的三维轴线，而对于图像序列中的每一后续时刻，则均以其前一时刻的3D血管轴线作为snake模型的初始位置，通过snake变形获得三维血管轴线，从而完成整个序列血管轴线的四维重建，具体步骤如下
a、以心电信号为同步信号，采集两个角度的、覆盖一个或多个心动周期的单面X射线冠状动脉造影同步图像序列，并对原始图像进行初步的滤波、去噪、畸变校正处理，增强视觉效果；
b、建立X射线血管造影系统在两个角度的投影模型，推导两幅不同角度的造影图像之间的几何变换关系
设点s1和s2为两次造影过程中X射线源焦点的位置，分别以s1和s2为原点，建立空间坐标系s1x1y1z1和s2x2y2z2；坐标系U1V1O1和U2V2O2为成像平面坐标系；D1和D2分别是s1和s2到各自成像平面的垂直距离，空间血管上的点P在成像平面A和B上的投影点分别为p1(u1，v1)和p2(u2，v2)，
则s1x1y1z1和s2x2y2z2之间的几何变换关系为
其中，R3×3为旋转矩阵R＝RY(α2)·RX(β2)·RX(-β1)·RY(α1)
为平移向量
L1和L2为X射线源焦点s1和s2到旋转中心的距离；(α1，β1)和(α2，β2)分别是A和B的造影角度；(x1，y1，z1)和(x2，y2，z2)分别表示空间点P在坐标系s1x1y1z1和s2x2y2z2中的坐标；
根据透视投影成像的几何关系，空间点的三维坐标与其二维投影坐标之间的关系为
其中，
c、在序列中第一时刻的图像中，通过人工取点获得感兴趣血管段在左右投影中的近似中心线，用折线表示，并利用外极约束得到两个角度间对应点的匹配；
d、对手动选取的采样点进行三维重建；
e、连接各3D点，所得折线作为3D snake的初始位置，通过使预先设定的能量函数最小，snake在空间中变形，最终得到具有最小能量的最优位置，就是第一时刻的3D血管轴线；
f、后续时刻血管轴线的四维重建
对于序列中的后续时刻，以前一时刻的3D血管轴线作为当前时刻snake的初始位置，通过snake变形，完成整个序列中血管轴线的四维重建。
2、根据权利要求1所述的冠状动脉血管轴线的四维重建方法，其特征是，所述X射线冠状动脉造影图像序列的两个采集角度之间夹角的取值范围为60°～120°。
3、根据权利要求1或2所述一种冠状动脉血管轴线的四维重建方法，其特征是，在第一时刻图像的感兴趣血管段中选取的三维重建采样点包括起点、终点和3～6个中间点。
全文摘要
一种冠状动脉血管轴线的四维重建方法，属医学检测技术领域。其技术方案是首先采集两个角度的X射线冠状动脉造影图像序列、建立X射线血管造影系统在两个角度的投影模型，并推导出两个角度图像之间的几何变换关系，再对手动选取的第一时刻图像中感兴趣血管段的采样点进行三维重建，并以连接各重建点所得折线作为初始位置，通过snake变形获得该段血管第一时刻的三维轴线，而对于图像序列中的每一后续时刻，则均以其前一时刻的3D血管轴线作为初始位置，通过snake变形获得三维血管轴线，从而完成整个序列血管轴线的四维重建。本发明大大减少了操作者参与所引入的不确定性和误差，提高了结果的可重复性，且操作简便、效率高。
文档编号G06T15/00GK101283911SQ200810055038
公开日2008年10月15日 申请日期2008年6月5日 优先权日2008年6月5日
发明者正 孙 申请人:华北电力大学