基于图像的回顾性脱机门控方法优化方案

考虑到ECG门控方法的缺点，可采用数字图像处理和分析技术，直接从连续回撤导管获取的ICUS图像序列中提取出隐含的心脏时相信息，然后根据该信息选择门控帧，即基于图像的回顾性脱机门控技术。该技术无需专用的门控图像采集装置，也无需利用心电信号，可解决没有同步记录ECG信号的ICUS图像序列的门控问题。

（1）方法概述

其示意图如图3-71所示，该技术包括以下步骤：首先，从图像序列中提取出反映心动周期的信号；然后，对该信号进行滤波，去除由非心脏运动因素引起的成分；最后，根据滤波之后的信号提取出特定心脏时相，并据此完成对门控帧的选择。选取门控帧的原则是：每个心动周期选择一帧；选择在每个心动周期（R-R间期）的舒张末期采集的图像。

对于第三个步骤，目前各方法均采用信号的极值作为各心动周期的采样点。因此各方法之间的主要差别在于前两个步骤，即隐含心动周期信息信号的提取，以及进一步的滤波方法。

针对第一个步骤，目前主要有两类方法下面具体介绍：

1）分析图像序列中血管结构或管腔尺寸的变化，具体如下：

例如，Nadkarni等^［70］首先对各帧图像进行手动分割，提取出管腔和血管壁的轮廓线，然后通过计算各帧图像中管腔轮廓线所包围区域的面积，分析血管腔横截面积的变化，从ICUS序列中选择出在近似相同相位处采集的各帧图像。该方法的应用前提是需要对各帧图像进行准确的分割，而由于目前还没有一种全自动的、鲁棒的血管内超声图像分割方法，因而限制了其应用。Zhu等^［71］通过计算各帧图像中一个手动设定的圆形感兴趣区域内的图像灰度变化，来分析管腔尺寸的变化。Rosales等^［69］根据管腔轮廓的椭圆近似，计算相邻帧图像中管腔轮廓之间的旋转角和重心之间的位移。Barajas等^［72］在由原始射频信号重建图像的过程中，通过分析管壁边界的变化，提取反映心脏运动的信号。此类方法的应用前提是需要对各帧图像进行准确的分割，分析结果的精度很大程度上取决于图像分割的精度。

图3-71 ICUS图像序列的回顾性脱机门控示意图

2）分析图像灰度的变化情况，具体如下：

通过分析图像灰度变化，完成对心脏时相的提取或门控帧的选择，无须对图像进行预先的分割、特征提取和跟踪等高层次的图像分析。此类方法的关键是选择适当的相似度（或差异度）测量标准，衡量不同图像之间的相似度或差异度。目前的方法都是采用涉及图像中所有像素的全局相似性测量。例如，Winter等^［73］提出的“Intelligate”方法降低了应用前提，即无须对ICUS图像进行预先的分割、特征提取和跟踪等高层次的图像分析，而是将两帧图像灰度特征的归一化互相关（Normalized Cross-Correlation，NCC）作为对ICUS序列中各帧图像之间相似度的衡量，然后将对门控帧的选择问题转化为对序列中各帧的归类，从而搜索出在各心动周期的舒张末期采集的图像。Malley等^［74］采用文献［73］的衡量各帧图像差异度的方法，构建整个图像序列的差异矩阵，并通过分析差异矩阵的周期性结构，估计出病人心动周期长度的近似值，再通过在差异矩阵中搜索具有最小累积差异值的路径，寻找在各心动周期的心脏运动最慢点处采集的一帧图像，组成门控序列。其不足之处在于，当病人的心率发生变化时可能得到错误的结果。文献［75］和［76］对［74］的方法进行了改进，前者为了构建更为精确的差异矩阵，首先对每帧ICUS图像进行特征提取，得到其图像描述子，然后采用图像描述子之间的距离作为对两帧图像差异度的衡量。这种方法虽然提高了差异矩阵的精度，但同时也提高了方法的复杂度。后者保证了算法在病人的心率发生变化时也能得到满意的结果。此类方法的主要不足在于采用全局相似性度量的分析结果极易受到图像中的噪声、纹理变化和背景区域的影响。

第二个步骤的关键在于滤波器的设计。文献［70］、［71］和［72］分别采用Butter-worth、Gabor和Daubechies滤波器来对反映心脏运动的信号进行滤波。文献［73］分别采用Butterworth和Gabor滤波器对信号进行滤波，并通过实验证明采用Gabor滤波器的滤波结果更满意。文献［77］探讨了带通滤波器的形状对采样精度的影响，并讨论了Gauss和Butter-worth滤波器的性能。

（2）方法举例

在连续回撤超声导管采集ICUS图像序列的过程中，伴随着周期性心脏运动的管腔形状变化会导致各帧管腔截面图像灰度特征的周期性变化，因此通过分析不同帧ICUS图像之间的灰度变化规律，可获得心脏的时相信息，并选取出在相邻心动周期中的相同相位处采集的图像，组成新的序列。基于上述思路，可按照如下步骤完成回顾性脱机门控：首先，分析各帧ICUS图像的灰度特征，构造差异矩阵；然后，从差异矩阵中找到累计差异值最小的路径，从而为各帧找到其在相邻心动周期中的对应帧；最后，选择在各心动周期的舒张末期采集的图像。各步骤的具体方法和对临床图像的实验结果举例如下：

1）构造ICUS序列的差异矩阵，具体如下：

对于匀速连续回撤超声导管采集的由n帧图像组成的ICUS图像序列｛I₁，I₂，...，I_n｝，首先分析各帧图像之间灰度特征的差异度，构造一个n×n维的差异矩阵D＝（d_i，j）（i，j＝1，2，..，n），其中d_i，j表示第i和j帧图像之间的差异值。选择两帧图像灰度值的归一化互相关作为对其差异度的衡量。设序列中各帧图像的大小均为M×N像素，对于I_i和I_j（i，j＝1，2，..，n），差异度d_i，j定义如下：

式中，μ_i和μ_j分别是I_i和I_j的平均灰度值。d_i，j越大，则I_i和I_j的差异越大。如此构成的差异矩阵具有如下特点：所有元素均为非负值，即d_i，j∈［0，1］；主对角线上的所有元素都是0，即i_i，i＝0；矩阵关于主对角线对称，即d_i，j＝d_j，i。

将D中的每个元素d_i，j用灰度值为255.d_i，j的像素点来表示可得到一幅n×n像素的灰度图像，其中像素越亮则表示对应的两帧图像的差异越大，反之则越小。图3-72是一个由248f组成的ICUS灰阶图像序列的差异矩阵，可以看出该矩阵具有明显的周期性结构，这是由于心脏的周期性运动所造成的ICUS图像内容的变化比其他因素引起的变化（包括血管本身的几何形态、搏动的血流引起的血管壁变形等）要快得多。

图3-72 包含248f图像的ICUS图像序列的差异矩阵

2）确定相邻心动周期中的对应帧，具体如下：

由于在相邻心动周期的近似相同相位处采集的图像应十分相似，它们之间的差异度也应呈现局部最小值。因此，通过在差异矩阵中寻找一条累计差异度最小的最优路径，可为ICUS序列中的各帧图像找到其在相邻心动周期中的对应帧。

为了找到满足要求的最优路径，首先须解决以下两个问题：确定最优路径的起点和确定搜索范围。第一个问题即确定序列的首帧I₁与第二个心动周期中的第几帧是对应帧。差异矩阵D中第一行的各元素｛d_1，1，d_1，2，d_1，3，...，d_1，n｝表示I₁与序列中各帧之间的差异度，且d_1，1＝0。曲线d_1，i～i，i＝1，2，...，n具有近似周期性的形状，各个谷点（局部极小值）所对应的i值间距近似相等，这是由心脏的周期性运动所造成的。随着i的增大，各局部极小值呈递增的趋势。也就是说，随着导管的回撤，在后续心动周期的相同相位处采集的图像与I₁的差别将越来越大。图3-73是对图3-72中的差异矩阵计算出的d_1，i～i曲线，显然第一个谷点（也是该曲线除d_1，1＝0之外的全局最小值）所对应的i值就是第二个心动周期中与I₁相对应的帧的序号，记为F。那么最优路径的起点即为d_1，F。由图3-73可知，F＝25f，表示该序列中第1和25f是在相邻心动周期的相同相位处采集的。

图3-73 由图3-72中的差异矩阵计算出的d_1，i～i曲线

第二个问题即确定相邻心动周期的对应帧之间，其序号之差的范围。由于在采集图像的过程中，是以匀速回撤超声导管，且图像的帧采样率也是固定的，因此在病人心率恒定的情况下，相邻心动周期中对应帧序号之差应等于以帧为单位的心动周期的长度。但是，考虑到病人心率可能发生变化的情形，可首先估计以帧为单位的心动周期长度的近似值C_C；然后，对于ICUS序列中的各帧，将对其在相邻心动周期中对应帧的搜索范围Δ设定为C_C／2≤Δ≤2C_C。

ICUS序列中间隔为if的两帧图像I_m和I_m＋i（m＝1，2，…，n-i）的平均差异度值D（i）如下：(www.xing528.com)

得到函数D（i）～i，i＝0，1，2，...，n-1，且 （0）＝0。图3-74a是对图3-72中的差异矩阵计算出的 （i）～i曲线，与曲线d_1，i～i类似，该曲线同样具有近似周期性的形状，且心率的近似值R（单位：次／min）应等于该曲线的重复频率。若图像采集速率为30f／s，则心动周期长度的近似值C_C＝（60×30）／R（f）。由图3-73b可知，R＝65（次／min），则C_C≈27（f）。

图3-74 由图3-72中的差异矩阵计算出的平均差异值曲线

a）时域波形（横轴为帧序号，纵轴为平均差异值） b）幅度谱（横轴是心率（单位：次／分钟），纵轴是频谱幅度）

采用动态规划（Dynamic Programming，DP）算法搜索差异矩阵中累计差异度最小的最优路径，不仅可保证解的全局最优性，而且具有计算效率高和数字计算稳定的优点。图3-75是n×n维的差异矩阵示意图，平面中的每个栅格点都表示一个数对［I_i，I_j］，对应两帧图像I_i和I_j，差异值为d_i，j。L₁是主对角线，其上所有点的值都是0。L₂是平行于L₁并且与其水平距离为Δ的直线。需要在该平面中寻找一条以d_1，F为起点、且经过的所有点的差异度累计值最小的路径。如前所述，将搜索范围限制在L₁和L₂之间宽度为Δ（C_C／2≤Δ≤2C_C）的条带区域内，即对于第r帧I_r，寻找其在下一个心动周期中的对应帧I_k时，须使r＜k≤r＋Δ。具体搜索步骤如下：

图3-75 n×n维差异矩阵示意图

在第r步（r≥3）中，确定终点为［I_r，I_k］（r＜k≤r＋Δ）、具有最小累计差异值的路径。方法是连接该点与前一步中的各点［I_r-1，I_I］（r-1＜I≤r-1＋Δ），计算出各条路径的累计差异值，并选择出具有最小累计差异值的路径。令c（r，k，r-1，I）为连接点［I_r，I_k］和［I_r-1，I_I］的路径的累计差异值，其公式如下：

c（r，k，r-1，I）＝d_r，k＋d_r-1，I （3-69）C（r，k）是终点为［I_r，I_k］的路径的最小累计差异值。将c（r，k，r-1，I）加在C（r-1，I）上，得到累计差异值Q（r，k，r-1，I）。对于步骤r中的每个节点［I_r，I_k］（r＜k≤r＋Δ），共有Δ个累计差异值。C（r，k）定义如下：

式中，C（2，I）＝d_2，I＋d_1，F；I＝3，4，…，2＋Δ。产生C（r，k）的节点［I_r-1，I_u（r-1）］成为［I_r，I_k］的前趋，这样属于终点为［I_r，I_k］的、具有最小累计差异值的路径的全部节点就是｛［I₁，I_F］，［I₂，I_u（2）］，...，［I_r，I_k］｝。重复上述过程，直到r＝n-Δ时得到的具有最小累计差异值的路径也就是整个问题的最优路径｛［I₁，I_F］，［I₂，I_u（2）］，...，［I_n-Δ，I_u（n-Δ）］｝，在图3-74中用I表示，它表示ICUS序列中，第1帧与第F帧、第2帧与第u（2）帧、…、第n-Δ帧与第u（n-Δ）帧是对应帧，其中u（.）是匹配函数。

为了减少计算量，还可对搜索过程附加约束条件，例如：唯一性约束，∀r，∃kI_r→I_k，保证只考虑各帧在相邻心动周期中的最佳对应帧；单调性约束，如果I_r→I_k并且I_r＋1→I_k′，那么k′＞k，限制可能解的搜索空间，避免“多对一”，即一帧与其相邻心动周期中的多帧相对应。

3）选择在舒张末期采集的图像，具体如下：

目前临床采用的ECG门控图像采集方法通常采用R波作为门控帧的标志，由于心脏在R波出现时处于一个心动周期中运动最缓慢的阶段，因此为了得到与ECG门控相似的结果序列，还需在最优路径I上找到表示在相邻心动周期的心脏运动最慢点处采集的图像对的点。若最优路径I经过d_i，j，则表示心脏在i和j两个时刻处于近似相同的状态，差异矩阵在点（i，j）呈现局部最小值；若第i和j帧同时又是在心脏运动最慢点采集的，那么在过点（i，j）、垂直于主对角线的方向上还应存在差异值很小的点。

首先采用45°高斯模板对差异矩阵进行滤波，强化具有上述特点的结构，得到滤波之后的差异矩阵D′。然后，在D′中，找到路径I上的最大值点（p₀，q₀）（在图3-76中用“★”表示），它表示第p₀和q₀帧是在相邻心动周期的心脏运动最慢点处采集的。以该点为起点，在I上找到表示其他对应帧的点，包括向下和向上搜索两个阶段，流程如图3-77所示，其中步骤2中的交点（q_i，j）和（j，p_i）在图3-76中都用“▲”表示。

图3-76 在滤波后的差异矩阵中搜索在舒张末期采集的各对应帧

临床采集的ICUS图像序列的门控结果如图3-78所示，从图中可以直观地看出，脱机门控序列较ECG门控结果更为平滑。

图3-77 搜索门控帧采样点的流程图

a）向下搜索 b）向上搜索

图3-78 ICUS序列的门控结果

a）原始ICUS序列纵向视图 b）基于图像的脱机门控结果 c）离线ECG门控结果

上述回顾性脱机门控方法不需要专门的心电门控图像采集装置，也无需利用ECG信号，同时亦无需操作者的手动参与，因而结果的客观性和可重复性高。由于采用了动态规划等高效算法，还可保证较低的运算量和较高的运算效率。