基于平均帧间差异度的门控方法优化

由于心脏的周期性运动，连续回撤导管采集的各帧IC-OCT图像的灰度特征也会发生周期性变化。由此，应首先构建图像序列的差异矩阵，并据此估算出平均心率的近似值；其次，提取出反映心动周期的一维信号，并对其进行带通滤波，去除由非心脏运动因素所致的成分；最后，通过检测滤波后信号的局部极值，提取出门控帧，完成对图像序列的回顾性脱机门控。下面介绍具体步骤。

（1）图像预处理

为了减小后续处理的数据量，提高算法效率，首先对IC-OCT图像进行灰度化处理，后续章节中的处理都是针对转换后的灰度图像进行的。

将彩色图像转化为灰度图像的方法主要有：单一分量法、最大值法、平均值法、加权平均值法、HLS模型算法、基于彩色空间距离的算法、基于梯度域的算法等。考虑到图像转换的保真度，可采用如下算法将RGB格式的彩色IC-OCT图像转换为灰阶为256的灰度图像，公式如下：

g＝R／2²＋R／2⁴-R／2⁶＋G／2＋G／2⁴＋G／2⁵＋B／2³-B／2⁶ （4-1）

式中，g为灰度图像中某个像素的灰度值；R、G、B分别为彩色图像中对应像素的红、绿、蓝分量。

图4-18是由300帧组成的图像序列灰度化后的结果，原始图像由美国LightLab公司生产的光学干涉断层成像系统采集，成像导管的回撤速率为1mm／s，帧速率是30f／s，图像的大小为256×256像素，灰度化以后的灰阶范围为［0，255］。其中图4-18a为纵向视图，可见血管壁呈现明显的锯齿状边缘；图4-18b为连续的8帧横向视图，显然不同帧之间的血管腔横截面存在平移和旋转。

图4-18 由300帧组成的图像序列灰度化后结果

a）纵向视图 b）横向视图

（2）估算平均心率

首先，对IC-OCT图像序列进行逐帧比较，利用两帧图像灰度值的归一化互相关计算差异值，建立差异矩阵，计算各帧图像之间的差异值。具体步骤如下为对于匀速连续回撤导管采集的由n帧图像组成的IC-OCT图像序列，计算两帧图像之间灰度值的归一化互相关（NCC），进而得到各帧图像之间的差异值d_ij如下：

式中，i，j＝1，2，…，n；F_i和F_j分别表示图像序列中的第i帧和第j帧，其大小均为p×q像素，平均灰度值分别为μ_i和μ_j。由式（4-2）可知，d_i，j∈［0，1］，且d_i，j越大，F_i和F_j的差异越大。由此构成的差异矩阵的特点为主对角线上所有元素的值都是0，即d_i，i＝0；所有元素的值均为非负值；矩阵关于主对角线对称，即d_i，j＝d_j，i。

然后，计算间隔为k帧的两帧图像F_m和F_m＋k（m＝1，2，…，n-k）的平均差异值D（k），公式如下：

得到曲线D（k）～k，其中k＝0，1，…，n-1，d_m，m＋k是第m帧和第m＋k帧图像的差异值， （k）∈［0，1］且 （0）＝0。由于心脏的周期性运动，该曲线具有近似周期的形状，其重复频率，即其幅度谱峰值所对应的频率值即为病人平均心率的近似值R（单位：次／min）。由于人的心率范围一般为45～200次／min，所以对 （k）的幅度谱曲线峰值的查找范围限定为［45，200］。

图4-18所示图像序列的差异矩阵如图4-19所示，矩阵中的每个元素用灰度值为255. d_ij的像素表示，得到一幅大小为300×300像素的灰度图像，其中像素越亮，表示对应的两帧图像差异越大，反之像素越暗，则两帧图像的差异越小。以帧间隔k为横轴，平均差异值为纵轴，得到平均差异变化曲线，如图4-20所示，该曲线的近似周期性形状能够明显反映出差异矩阵的周期性变化。平均差异曲线的幅度谱曲线如图4-21所示，在［45，200］范围内有一明显峰值，此峰值所对应的频率即近似为平均心动速率。

图4-19 300帧IV-OCT图像的差异矩阵

图4-20 300帧IV-OCT图像的平均差异曲线

图4-21 平均差异值曲线的幅度谱曲线

图4-22 一帧IC-OCT图像(www.xing528.com)

（3）提取隐含心动周期的信号

一帧IC-OCT图像中背景区域内像素的灰度值在整个图像序列中基本不发生变化，所以这些像素不包含心脏运动信息，为方便起见，将其称为静止像素。为减小后续步骤的计算量，需先将它们去除。

除静止像素外，血管区域内的像素分为两类：一类为包含心脏运动信息的像素，大部分分布于血管壁上，例如图4-22中的黄色十字所示的像素，称为动态像素；另一类同时包含了心脏运动信息和非心脏运动信息（如呼吸和身体的移动等），称为干扰像素，在提取心动周期信号时，这类像素也需要被滤除。具体步骤如下为首先，对于一幅IC-OCT图像中的每个像素，统计序列中各帧图像上相应位置处的灰度值，产生一个一维灰度变化信号g_i，j（m），其中m表示帧序号，i，j＝0，1，…，N-1分别表示像素的横、纵坐标。如果一帧IC-OCT图像的大小为N×N像素，那么共产生N²个一维灰度变化信号。分析每个一维灰度变化信号的频谱，对于静止像素，由于其灰度值自始至终基本没有变化，所以其一维灰度变化信号近似为直流信号，其幅度谱在平均心率点处没有明显的波峰；对于被干扰的血管区域像素，其灰度值的变化没有明显的规律，所以其幅度谱曲线在心率点处的波峰也不明显，其他频率分量多且幅度较大；对于动态像素，其幅度谱曲线在平均心率点处有明显的峰值，其他频率分量的幅度明显小于波峰。

然后，提取由动态像素产生的一维灰度变化信号。设定判决条件为在平均心率处的频谱幅值大于其他频率处幅值的三倍，保留符合判决条件的一维灰度变化信号，舍去不符合的。最后剩余信号的频谱取平均，并作逆傅里叶变换，即可得到一个隐含心动周期的一维信号g（m），其中m表示帧序号。

所有像素平均后（即提取动态像素前）的一维灰度变化信号的时域及幅度谱曲线如图4-23所示，对动态像素的一维灰度变化信号求平均后的时域及幅度谱曲线如图4-24所示。由于提取后的信号在计算时已滤除直流分量，所以图4-24a的灰度值为相对值，不具有参考意义。通过对比幅度谱可以看出，提取后动态像素得到的一维信号频谱在平均心率点处有明显的峰值，其他频率处的幅值均远小于该峰值，去除了干扰像素的影响。

图4-23 提取动态像素前的一维信号

a）时域波形 b）幅度谱曲线

（4）隐含心动周期信号的滤波

隐含心动周期的信号g（m）中仍包含很多非心脏运动因素，如：管腔本身不规则的几何形状、不稳定的心率以及呼吸运动等，在频谱上表现为除心率外的其他频率分量。本书应用带通滤波器来滤除由非心脏运动因素所致的频率分量，最终得到信号g^（f）（m）。

带通滤波器的参数设计如下：中心频率为前述步骤中估算出的平均心率R；通带范围为［R-0.3R，R＋0.3R］，这是由于在图像采集过程中，病人的心率会以R为中心上下波动，上述范围既能使心率变化的范围包含在通带内，又能滤除非心脏运动因素的干扰。

图4-24所示的一维信号进行滤波后的时域波形与幅度谱曲线，如图4-25所示。由图像序列估算出的平均心率为72.24Hz，因此设置带通滤波器的通带上限截止频率为50.59Hz，下限截止频率为93.91Hz。可以看出，250Hz以上的高频分量基本全部滤除，心率的二次谐波幅度降低为心率幅度的近十分之一，表明非心脏运动分量基本全部滤除。在时域波形中，不同局部极值的平均灰度值不同，说明在滤除噪声的情况下，保留了不同时刻心动周期的差异性。

图4-24 提取动态像素后所得的一维信号

a）时域波形 b）幅度谱曲线

图4-25 滤波后的一维信号

a）时域曲线 b）幅度谱曲线

（5）选择门控帧

在不同心动周期的同一相位处采集图像，可以得到最小的帧间差异值。由于信号g^（f）（m）的波峰或波谷相对于其他位置更容易识别，因此对信号g^（f）（m）进行极值检测，提取出门控帧序列。

图4-18所示图像序列的门控结果如图4-26所示。为了便于比较，将门控后纵向视图人为拉伸至原序列长度。通过门控前后的对比可以看出，门控后纵向视图的视觉效果得到了很大的改善，抑制了血管壁的锯齿效应，管壁变得平滑连续。

图4-26 IC-OCT图像序列的门控结果

a）门控前 b）门控后