CAG图像中的运动跟踪技术优化

早期的图像跟踪多数是采用模板匹配方法，通过对包含目标物体的窗口在图像中进行匹配实现对目标的跟踪。由于血管是狭长的结构，不管采用矩形窗口还是圆形窗口，血管都只占据窗口中很小的一部分，因此很难设定一个包含整个血管段的窗口。另一方面，如果选择一系列小窗口对血管段进行跟踪，由于窗口之间缺乏足够的约束，也很难得到满意的结果。

从20世纪80年代后期开始，变形模型在跟踪变形体的运动中开始广泛采用。例如，Kass等^［57］在1988年提出的活动轮廓模型（Active Contour Model，即snake），可以在不需要更多的先验知识或高层处理结果指导的情况下，实现自动追迹，提取动态图像中欲处理目标的闭合、光滑、连续的轮廓，并且具有较强的抗噪声能力。snake模型是一条可变形的参数曲线，对该曲线构造合适的变形能，并通过最小化该函数，使得曲线在内部力和外部约束力的作用下变形，最终精确收敛到目标的边缘，得到具有最小能量的目标轮廓。模型的能量函数采用积分运算，具有较好的抗噪性，对目标的局部模糊也不敏感。自出现以来，在原有模型的基础上又出现了许多改进的模型，已经被成功地应用于边缘提取、图像分割、运动跟踪等许多领域。

本节介绍应用活动轮廓模型对CAG图像序列中的冠脉血管骨架进行跟踪的方法，将序列中前一时刻的血管骨架作为下一时刻snake的初始位置，变形得到下一时刻的骨架。

（1）snake模型原理

snake模型实质上是一个能量最小化模型，snake是定义在待分割图像上的一条任意闭合曲线，对该曲线构造合适的变形能。其形变过程就是snake在外部能量（外力）和内部能量（内力）的作用下向目标轮廓靠近的过程。外力的吸引使snake不停向着真实轮廓移动，而内力在保持snake拓扑性的同时随着snake的移动而变化，最终达到内外力之和等于零（对应于能量最小）。这时，snake就停留在物体的真实轮廓上。这种动态逼近方法所求得的边界曲线具有封闭、光滑等优点。这个最小化过程一般由求解Euler-Lagrangian方程，用迭代算法得到。

经典的snake模型是一条参数曲线c（s）＝（x（s），y（s）），s∈［0，1］。模型的能量函数定义如下：

其中内部能量E_int定义如下：

式中，α（s）和β（s）为权重，其作用分别是调节snake的张力和刚性；c′（s）和c″（s）分别为c（s）的一阶和二阶导数。内部能量使得参数曲线保持平滑和收缩的性质，其中的一阶导数项代表参数曲线的弹性能量，控制轮廓的“应力”，使得曲线保持收缩，使曲线具有类似弹性绳的特性，在点数不变的情况下减小该项意味着参数曲线的不断收缩；二阶导数项是曲率，使曲线保持平滑，曲线具有一定的刚度，减小该项可以控制参数曲线趋于平滑。

外部能量函数E_ext是定义在整个图像平面I（x，y）上的标量函数，是保证活动轮廓收敛的外部力，反映了图像的某些特性，如边缘等。它决定着轮廓的移动方向，由于这些方向只能根据特定的问题而定义，因此E_ext没有统一的数学表达式，必须从问题本身的特征出发，利用灵活的图像处理技术得到。在Kass的snake模型中，外部能量由两部分组成，公式如下：

E_ext＝E_image＋E_con （2-86）式中，E_con表示外部限制作用力产生的能量，叫约束力；E_image表示图像作用力产生的能量，也叫图像力，是直线、边缘和界限能量的线性组合，其公式如下：

E_image＝W_lineE_line＋W_edgeE_edge＋W_termE_term （2-87）式中各项均可以从图像I（x，y）计算得到，W_line、W_edge和W_term为特征系数，E_line是线能量，等于图像亮度（灰度），E_line＝I（x，y），通过线性特征系数W_line的正负，可以便snake靠近亮线或暗线。E_edge是边缘能量，也可由一个非常简单的函数表达为E_edge＝|▽-I（x，y）|²，边缘特征系数W_edge控制对轮廓所在区域的强度梯度的约束，E_term是用高斯函数平滑过处理的图像中各级轮廓线的曲率，由特征系数W_term决定其影响。

求解snake模型就是要使E最小化。将式（2-84）取极小值，满足Euler方程的具有最小能量的一阶偏微分方程为

αc″（s）-βc″″（s）-E_ext＝0 （2-88）

Kass采用变分的迭代方法求解该方程。

snake模型具有统一的、开放式的描述形式，为图像中变形体跟踪技术的研究和创新提供了理想的框架。在实现活动轮廓模型时，可以灵活地选择约束力、初始轮廓、作用域等，以得到更佳的分割效果。所以snake模型自其问世以来，一直是图像变形体跟踪研究中的一个热点。

（2）首帧图像中血管段骨架的提取

应用snake模型必须要解决两个问题：一是要选择合适的能量函数；二是要确定适当的初始轮廓位置，否则snake将收敛于错误的结果。在进行图像序列的跟踪时，对于序列中的首帧图像，传统的方法是通过人工取点提取初始血管中心线，作为snake的初始位置，然后利用snake模型方法得到第一幅图像的真实的中心线。对于其后的序列图像，根据血管形状不突变的特点，把第n幅图像中snake的停留位置作为第n＋1幅图像中snake的初始位置，这样就能够实现序列图像的自动分割和跟踪，流程图如图2-73所示。

在目标轮廓上取n个点，snake能量函数的离散形式如下：

其中E int由式（2-85）确定，导数项可用离散点的差分近似表示。设活动轮廓上的各点为v_i（i＝1，2，…，n），α和β取常数值，采用有限差分方法有

图2-73 采用snake模型进行CAG图像中血管骨架的二维运动跟踪流程图

那么E_int可以表示为

式中的第一项用于限制snake上点的距离，使两点不能相距过远（或过近），相当于张力，使其趋于最小则可保持曲线的连续；第二项是一点与相邻两点所成夹角的函数，即曲率，使其趋于最小则可保持曲线的平滑。连续项约束 实质是snake上点间距离的缩小，在实际求解过程中会使曲线收缩，甚至收缩为图像中的一点。为此，可使用点之间的平均距离 与当前所处理的两点之间距离之差作为约束，即 ，使点平均分布，与平均距离接近的点将具有较小的能量值。这样既能满足连续性的要求，又不会产生snake收缩的现象^［58］。在每次迭代结束时， 的值被更新。E_int的第二项为对曲率的约束，即平滑性约束。由于连续项约束的作用是使snake上的点按等距离分布，因此， 乘以一个常数，即可得到对曲率的合理估计。综上所述，得到如下公式：

其中

(www.xing528.com)

snake点间距与总点数是相互联系的两个参数。对于一个特定目标，点间距增大，则点数减少；点间距减小，则点数增多。事实上，点间距的大小取决于处理问题的性质。对于造影图像中血管的提取问题，当血管的轮廓比较简单时，可以选取较大的间距；而当目标的轮廓比较复杂时，为了能以较高的精度对目标进行拟合，则需要较小的间距。当间距较小时，snake上的点数较多，则计算量增大。

在提取序列中第一幅图像的血管段中心线时，外部能量只包含图像力E_image分量，它由图像点的灰度和灰度梯度的组合构成，公式如下：

其中

式中，I_x和I_y分别是图像的灰度函数I（x，y）对横纵坐标的偏导数。对I（x，y）的每个像素，通过计算2×2邻域矩阵的平均有限差分来作为x和y方向的偏微分的一阶近似，公式如下：

由于血管造影过程中，当造影剂的浓度达到峰值时，所拍摄的图像上血管的灰度值应比背景的灰度值小，所以这里权重系数γ取正值。而且由于提取的是血管中心线，不是边缘，所以λ也取正值。

通过迭代运算，使式（2-89）表达的snake能量函数最小，snake在内力和外力的牵引下不断变形，最后停留在最优位置。对这些离散点进行三次B样条曲线拟合，得到一条连续曲线，采用较少的参数实现对血管中心线的子像素级精度描述。然后在曲线上再进行均匀采样，用这些采样点作为下一时刻snake的初始位置。

图2-74和图2-75是对两幅左冠CAG图像的血管段骨架跟踪结果，手动勾画出感兴趣血管段的初始位置，应用snake算法，使能量函数最小，迭代计算出该血管段中心线的真实位置。

图2-74 RAO角度的左冠造影图像及血管段骨架提取结果

a）原始图像 b）snake模型的初始位置 c）血管段骨架提取结果

（3）整个图像序列中血管段骨架的运动跟踪

在提取出序列中第一幅图像的血管中心线之后，对于其后的序列图像，根据血管形状不突变的特点，把第n幅图像中snake的停留位置作为第n＋1幅图像中snake的初始位置，这样就能够实现序列图像的分割和跟踪。

在对整个序列进行运动跟踪时，外部能量函数除了式（2-94）表示的图像力之外，还可增加相似度匹配，通过对相邻帧之间的目标进行配准实现对运动物体的跟踪。假设在足够短的时间间隔内，运动前后血管投影的灰度不发生变化，其公式如下：

图2-75 LAO角度的左冠造影图像及血管段骨架提取结果

a）snake模型的初始位置 b）血管段骨架提取结果

I（x＋Δx，y＋Δy，t＋Δt）＝I（x，y，t） （2-97）那么得到序列中第一帧图像的血管中心线之后，在开始对序列中其他帧进行运动跟踪时，外部能量函数可采用如下的形式：

式中，u_i是当前snake上的第i个点；p_i是前一帧图像中心线上的第i个点。这样，外部能量中除了包含图像力，即吸引曲线向着图像中灰度最小、灰度梯度变化最小的区域内移动之外，还包含外部约束，使得运动前后血管上对应点的灰度变化最小。

由于血管是开轮廓，因此应用snake进行血管跟踪，如果不对血管中心线的两个端点附加其他的外力约束，可能出现中心线沿血管收缩、伸长和漂移的现象。为了防止出现这些现象，操作者在选取血管段的时候，需选取两个特征点（即分叉点和端点）之间的血管段，那么在确定了血管段的起点和终点之后，为了避免snake收缩／伸长，需在外部能量中附加弹簧力，公式如下：

式中，K_spring是弹性系数；v₀和v_n-1分别是血管段的起点和终点。snake在变形过程中，其两端点是固定的，弹簧力约束使得与端点相邻的点偏离相应端点的距离不致太大，结合连续性和平滑性约束，最终得到光滑、连续的中心线。

对一个左冠CAG图像序列进行血管段骨架运动跟踪的实验结果如图2-76所示，其中图2-76a是首帧原始图像，选择出旋支上感兴趣血管段的起点和终点，连接两点的线段作为snake模型的初始位置。通过使能量函数最小，迭代计算出该血管段中心线的真实位置，并对离散的snake点进行三次B样条曲线拟合，结果如图2-76b所示。从图像序列中任选了10个时刻的图像，对感兴趣血管段的骨架进行运动跟踪，结果如图2-76c～图2-76e所示，从图中可直观地了解该血管段在整个心动周期中的形状和位置变化。由图2-76e可见，从时刻3～10是一个心动周期，在周期结束时，血管段会回到其在周期开始时的位置。

图2-76 对一个左冠CAG图像序列进行血管段骨架运动跟踪的实验结果

a）首帧图像中snake模型的初始位置 b）血管段骨架提取结果 c）第二帧图像中snake模型的初始位置 d）第二帧图像中血管段骨架的提取结果 e）图像序列中10个时刻的血管段骨架