一种抑制IVUS图像序列中运动伪像的方法

孙 正

（华北电力大学电子与通信工程系，河北 保定 071003）

血管内超声(intravascular ultrasound，IVUS)成像是临床诊断冠心病的一种基于导管的腔内成像技术，具有独特的在活体中观察血管壁、管腔及斑块形态性质，甚至管壁功能状态的特点，已经成为 X射线冠状动脉造影(X-ray coronary angiography)的弥补影像技术。

由于心脏运动以及血压变化所造成的运动伪像是影响IVUS图像的视觉效果，进而影响定量分析和血管三维重建精度的主要因素。目前主要采用门控的方法减少 IVUS序列中的运动伪像，包括3种方式：

一是心电(electrocardiogram，ECG)门控的图像采集方法：即采用 ECG门控的马达回撤超声导管，根据 ECG信号只在每个心动周期的相同时相（一般是R波）采集图像，改善心脏搏动所致的伪像[1]。该方法的缺点在于需要专门的ECG门控采集装置，而目前临床采用的多数IVUS系统不包含此功能。并且由于每个心动周期只采集一帧，因此与连续回撤导管的非门控采集方式相比，图像采集时间延长至少3倍[2]。

二是离线 ECG门控方法：即连续回撤超声导管采集图像序列，同时记录 ECG信号。待介入检查过程结束后，由医生根据 ECG信号从图像序列中选取在相同心脏时相采集的图像。这种方式虽然不会延长介入检查的时间，但是其不足之处在于：ECG信号是心脏的生理电特性通过骨骼、肌肉和皮肤等反映在体表的信号，而 IVUS图像是通过插入血管腔内的超声探头采集的，两者不可能完全同步；门控结果的客观性和可重复性差，受操作者的临床经验和主观因素影响较大；虽然多采用 0%点即舒张末期作为采样点，但无法确定采用该点是否可以获得最稳定的门控图像序列；目前尚无公认、明确的标准，来确定采用心动周期中的哪一时相作为最佳采样点，可获得最大帧间稳定度[3]。

三是基于图像的离线门控方法：基于 ECG门控技术在抑制运动伪像方面的不足，近年来，运用数字图像处理技术，从连续回撤超声导管采集的IVUS图像序列中提取出隐含的心动周期信息，并据此从每个心动周期中选取一帧图像，组成同步采样序列，即基于图像的门控(image-based gating)技术[2-8]，逐渐成为该领域的研究热点。这种方法虽然可以避免离线 ECG门控方法的缺点，但是由于每个心动周期只选择一帧，需要抛弃大量帧，因而可能会丢失很多有诊断价值的信息。

目前，对非门控IVUS图像序列直接进行运动估计的方法包括光流法(optical flow，OF)[9]和块匹配法(block matching，BM)[10]。此类方法既不需要门控图像采集装置，无需记录ECG信号，也不需要抛弃有用帧，可保证图像数据集合的完整性。但是光流法的估计结果极易受到血流随机运动的干扰，且对噪声十分敏感；块匹配法需要明显的标志物（例如钙化，血管分叉等），而此类标志并非在每帧中都存在。

本文在综合分析IVUS图像序列中由周期性心脏运动所致运动伪像的产生机制和表现形式的基础上，对覆盖多个心动周期的非门控 IVUS图像序列直接进行刚性运动伪像的定量估计和补偿。首先建立横向视图中的血管壁在心动周期中的刚性运动模型；然后定量估计相邻帧间的血管横断面刚性运动参数（包括位移和旋转角），并采用谱分析的方法分离运动分量与几何分量；最后通过补偿运动分量，达到抑制运动伪像的目的。无需利用 ECG信号，也不需要抛弃任何一帧图像，在抑制运动伪像的同时保证图像数据集合的完整性。

1 方 法

1.1 运动伪像的产生机制和表现形式

IVUS图像序列中由周期性心脏运动和搏动血流所致的运动伪像的主要表现形式为：

1） 帧间的错位，即横断面图像序列（图1(a)）中相邻帧之间血管横截面的平移和旋转（图1(b)）；

2） 在沿管腔长轴方向的纵向视图中，血管壁边缘呈现锯齿形（图1(c)）。

图2说明了运动伪像的产生机制。假设对一段均匀直径的直血管段进行IVUS成像，在回撤超声导管的过程中，导管始终位于各帧断层图像的中心。若无运动伪像，且假设导管回撤路径与管腔轴线重合，则将各帧图像按照采集顺序依次排列得到的纵向视图中，血管壁的上下轮廓应近似为直线（图 2(a)）。若血管处于周期性运动状态，则可能导致导管和管腔之间的相对运动。因而在回撤导管的过程中，采集到的各帧断层图像记录的并非血管的真实状态，在纵向视图中将出现“锯齿样”的血管壁上下轮廓，如图2(b)所示。这种伪像不仅会影响对血管腔和斑块形态的观察，同时也会给形态参数的准确测量和血管的三 维重建带来很大困难。

图1 IVUS图像序列

图2 IVUS图像序列示意图

1.2 建立血管壁刚性运动模型

IVUS断层图像序列中相邻切片之间管腔横截面重心的位移和空间方向的改变主要由运动和几何两方面的因素造成[3,11-13]：

1） 运动：即外部因素，指由周期性心脏运动和搏动的血流造成的超声导管相对于管腔的运动和血管形态的变化；

2） 几何：即内部因素，指血管腔本身不规则的几何形状。

从各帧断层图像中分割出血管腔边界之后，相邻帧之间管腔横截面的平移和旋转可分别用管腔边界曲线重心的位移(,)xyΔΔ和边界曲线之间的旋转角αΔ来表示，它们分别由两部分组成

其中，脚标d和g分别表示运动分量和几何分量。

如图3所示，设相邻两帧中的管腔边界曲线分别为1γ和2γ，其重心分别为，在运动和几何两方面因素的共同作用下，1γ和2γ之间的位移为，1γ和2γ之间的旋转角为Δα，旋转中心为导管(x0,y0)，即

其中，Δx、Δy和Δα如式(1)所示分别由运动和几何分量两部分组成。对于图像序列中的各相邻帧，本文通过计算出位移(Δx,Δy)和旋转角Δα，并从中分离出Δxd、Δyd和Δαd，完成对运动伪像的补偿。

图3 相邻帧IVUS图像运动参数示意图

1.3 运动伪像的估计和补偿

1.3.1 提取IVUS图像中的血管壁边界

对IVUS图像进行分割，提取血管内腔边界和血管壁中-外膜边界，不仅是本文抑制运动伪像方法的重要步骤，而且也是量化分析[14]、三维重建[15]以及 IVUS与其它冠状动脉图像融合[16]的基础和保证其精度的关键。近年来，提出了很多计算机辅助的自动分割方法[17-18]。考虑到处理效率等因素，本文采用一种三维分割方法[19]，与二维分割方法是对图像序列进行逐帧串行处理相比，三维分割技术可实现对图像序列的并行处理，提高处理效率。

1.3.2 计算运动参数

对于图像序列中的各相邻帧，例如第k-1和k帧（k=2,3,…,M，M是IVUS图像序列的总帧数），计算血管内腔边界曲线的几何中心作为对其重心的近似，并如§1.2中所述方法计算出重心之间的位移量和旋转角和kαΔ包含运动分量和几何分量)两部分，需从中分离出运动分量,dkxΔ、,dkyΔ和,dkαΔ。

1.3.3 分离运动分量和几何分量

本文对计算出的kxΔ、kyΔ和kαΔ（k=2,3,…,M）分别进行傅立叶变换，得到其幅度谱ΔXk、ΔYk和ΔAk。由于由血管本身的不规则几何形状所引起的相邻帧之间管腔横截面空间方向和重心位置的变化速度，远小于由周期性心脏运动所致的变化速度，因此在ΔXk、ΔYk和ΔAk中，高频分量对应于刚性运动分量，且其变化频率应等于心率，而低频分量对应于几何分量 ( Δxk,g, Δyk,g, Δαk,g)。分别对ΔXk、ΔYk和ΔAk进行高通滤波，则滤波器的输出即是刚性运动分量的幅度谱(ΔXk,d, ΔYk,d, Δ Ak,d)，对其进行逆傅立叶变换即可得到 ( Δxk,d, Δyk,d, Δαk,d)的估计值。

由于运动分量 ( Δxk,d,Δyk,d, Δαk,d)主要由周期性心脏运动引起，因此本文将高通滤波器的通带截止频率设定为病人的心率值（单位：次/秒，即Hz）。采用文献[6]中的方法从IVUS图像序列中估计平均心率的近似值，原理如下：

首先，对IVUS图像序列进行逐帧比较，通过计算两帧图像之间灰度值的归一化互相关(NCC)

计算各帧图像之间的差异值

式中，i,j= 1 ,2,… ,M；Ii和Ij分别表示IVUS图像序列中的第i帧和第j帧，其大小均为N1×N2像素，平均灰度值分别为μi和μj；M是图像序列的总帧数。

然后，计算间隔为k帧的两帧图像Im和Im+k（m= 1,2,… ,M-k）的平均差异值

得到曲线～k，k=0,1,2,… ,M-1，且= 0 。该曲线具有近似周期的形状，这是由心脏的周期性运动所造成的，该曲线的重复频率，即其幅度谱曲线峰值所对应的频率即为病人心率的近似值R（单位：次/分钟，即Hz）。

1.3.4 补偿运动伪像

2 实验结果与分析

为了验证方法的可行性并估计精度，本文分别利用计算机模拟的血管腔轮廓序列和临床采集的IVUS图像序列进行实验，并对实验结果进行分析。

2.1 模拟图像实验

图4 7帧模拟管腔轮廓序列

为了评价§1.3中所述对血管腔轮廓序列刚性运动参数的估计精度，本文将用封闭曲线表示的管腔轮廓序列按照已知的位移和旋转角进行平移和旋转，模拟生成运动参数已知的管腔轮廓序列。图4(a)和图4 (b)分别为7帧模拟管腔轮廓的横向视图和完成重心计算的纵向排列图，图4(c) ～ (e)为各偏移量的柱状图。

图5和图6分别为对20帧和35帧模拟管腔轮廓曲线序列进行刚性运动补偿前后的横截面视图和纵向视图。可以看出，完成运动补偿后的横截面视图中相邻帧之间的扭动大大减小，在纵向视图中体现为管腔边界的锯齿形尖峰与水平基线之间的距离大大缩小。

图6 35帧模拟图像实验结果

2.2 临床图像实验

临床图像序列是采用Jomed Endosonic超声成像仪采集的，采用2.9F 30MHz单轨机械超声探头，由马达驱动自动回撤探头导管，回撤速度为0.5mm/s，当探头以1800转/分作360°旋转时，可以30帧/秒的速率连续获得血管横轴实时切面图像。

图7 62帧临床图像实验结果

图7和图8分别为62帧和80帧临床图像序列的运动伪像补偿前后的横向视图（血管腔边界曲线）及纵向视图。可以看出，完成运动伪像抑制后的横截面视图中相邻帧管腔轮廓之间的扭动大大减小，在纵向视图中体现为外膜轮廓的锯齿形边缘变的相对平滑，因而很大程度上抑制了图像序列中存在的运动伪像。

图8 80帧临床图像实验结果

由于目前还没有一种客观、定量评价抑制IVUS运动伪像方法精度的体系[6]，因此，我们将采用本文方法进行运动伪像抑制的结果和采用基于图像的离线门控[8]结果进行比较，评价本文方法的精度。图9是分别对120帧和200帧的临床图像序列进行实验的结果，可见虽然基于图像的离线门控法可实现对运动伪像的抑制，但是门控序列的长度相对于原序列而言明显变短，即丢弃了大量的帧数据。而本文方法在完成运动伪像抑制的同时，不丢弃任何一帧数据，保证了原始图像数据集合的完整性。

2.3 误差源分析

本文的运动伪像抑制方法包括图像的采集和预处理、血管腔轮廓提取、管腔轮廓重心及偏移量的计算、运动分量的滤出、对运动分量的补偿等步骤。除去计算误差之外，误差主要来自两方面：

一是图像的分割：由于本文方法是直接对分割出的管腔轮廓进行运动估计的，因而IVUS图像的分割精度直接影响对运动伪像的补偿效果。因此，进一步提高二维轮廓提取的准确性有助于提高运动伪像的补偿精度。

二是滤波器的设计：由于本文采用高通滤波的方法分离由心脏的周期性运动引起的偏移量，因而高通滤波器的设计（包括原型的选择和参数的设定）直接影响对运动伪像的抑制效果。考虑到运算复杂度等因素，本文选择了 Butterworth滤波器，并将通带截止频率设定为从IVUS图像序列中估计出来的平均心率值。对于滤波器原型的选择，文献[3]首次探讨了滤波器的通带形状对采样精度的影响，并提出采用统计的方法确定哪种滤波器更适合于心脏时相信息的提取。在后续工作中，我们将着重分析滤波器的性能对于运动伪像抑制效果的影响。对于平均心率值的估算，若图像采集过程中同时记录了 ECG信号，则可直接利用ECG记录的心率值；若无ECG信号可供利用，则采用文献[6]的方法，在病人心率保持恒定的情况下，可得到精确的估算结果。但是若病人的心率发生波动，则其估计结果存在一定误差，详细讨论参见文献[6]。

图9 临床IVUS图像序列运动伪像抑制实验结果

3 结 论

本文提出了一种抑制冠状动脉内超声图像序列中由心脏运动所致运动伪像的方法，与现有的心电门控方法（包括联机和脱机两种方式）相比较，既不需要专门的心电门控图像采集装置，也不需要同步记录心电信号，而是运用图像分析技术，对连续回撤超声导管采集到的IVUS图像序列数据进行处理，因而适用于没有记录 ECG信号的IVUS记录序列。与基于图像的门控方法相比，本文方法不需要抛弃图像序列中的有用帧，可在抑制运动伪像的同时，保证图像数据集合的完整性。

[1]胡广禄, 陈兴新. 心电门控技术及其应用[J]. 医疗装备, 2000, 13(1): 17.

[2]de Winter S A, Hamers R, Degertekin M, et al. A novel retrospective gating method for intracoronary ultrasound images based on image properties [C]//Proceedings of IEEE International Conference on Computers in Cardiology, Kassandra Chalkidiki,Greece: IEEE Press, 2003, 30: 13-16.

[3]Hernàndez-Sabaté A, Gil D, Garcia-Barnés J, et al.Image-based cardiac phase retrieval in intravascular ultrasound sequences [J]. IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control,2011, 58(1): 60-72.

[4]Zhu H, Oakeson K D, Friedman M H. Retrieval of cardiac phase from IVUS sequences [C]//Proceedings of SPIE Conference on Medical Imaging 2003:Ultrasonic Imaging and Signal Processing, Bellingham WA: SPIE Press, 2003, 5035: 135-146.

[5]Nadkarni S K, Boughner D R, Fenster A. Image-based cardiac gating for three-dimensional intravascular ultrasound [J]. Ultrasound in Medicine and Biology,2005, 31(1): 53- 63.

[6]O’Malley S M, Granada J F, Carlier S, et al.Image-based gating of intravascular ultrasound pullback sequences [J]. IEEE Transactions on Information Technology in Biomedicine, 2008, 12(3):299-306.

[7]Carlo G, Oriol P, Oriol R L, et al. Robust image-based IVUS pullbacks gating [C]//Proceedings of 11th International Conference on Medical Image Computing and Computer-Assisted Intervention Part II(MICCAI 2008), New York City, USA: MICCAI Society Press, 2008, 5242: 518-525.

[8]孙 正. 抑制冠状动脉内超声图像序列运动伪影的离线门控方法[J]. 光电子·激光. 2010, 21(4):632-638.

[9]Danilouchkine M G, Mastik F, Steen A F W. Accuracy in prediction of catheter rotation in IVUS with feature-based optical flow——a phantom study [J].IEEE Transactions on Information Technology in Biomedicine, 2008, 12(3): 356-365.

[10]Danilouchkine M G, Mastik F, Steen A F W. Two methods for catheter motion correction in IVUS palpography [C]// Proceedings of 2007 IEEE International Ultrasonic Symposium (IUS), New York NY: IEEE Press, 2007: 1724-1727.

[11]Nadkarni S K, Mills G, Boughner D R, et al. A pulsating coronary vessel phantom for two-and threedimensional intravascular ultrasound studies [C]//Proceedings of the 4th International Conference on Medical Image Computing and Computer-Assisted Intervention (MICCAI 2001), Utrecht, Netherlands:MICCAI Society Press, 2001, 2208: 1164-1165.

[12]Holzapfel G, Gasser T, Stadler M. A structural model for the viscoelastic behavior of arterial walls:continuum formulation and finite element analysis [J].European Journal of Mechanics A-Solids, 2002, 21:441-463.

[13]Rosales M, Radeva P, Rodriguez-Leor O, et al.Modeling of image-catheter motion for 3-D IVUS [J].Medical Image Analysis, 2009, 13(1): 91-104.

[14]Baldewsing R A, Danilouchkine M G, Mastik F, et al.An inverse method for imaging the local elasticity of atherosclerotic coronary plaques [J]. IEEE Transactions on Information Technology in Biomedicine, 2008, 12(3): 277-289.

[15]Paul S. 3-D Intravascular ultrasound (IVUS) and IVUS-Palpography: insights into the mechanical behavior of the coronary vessel wall [J]. The International Journal of Cardiovascular Imaging,2006, 22: 153-155.

[16]Sun Zheng. Three-dimensional reconstruction of vessels from intravascular ultrasound sequence and X-ray angiograms [C]//Proceedings of 2nd International Conference on BioMedical Engineering and Informatics (BMEI 2009), Tianjin,China: IEEE Press, 2009, 1: 444-448.

[17]Bovenkamp E G P, Dijkstra J, Bosch J G,et al.User-agent cooperation in multiagent IVUS image segmentation [J]. IEEE Transactions on Medical Imaging, 2009, 28(1): 94-105.

[18]Zhang Qi, Wang Yuanyuan, Wang Weiqi, et al.Automatic segmentation of calcifications in intravascular ultrasound images using snakes and the contourlet transform [J]. Ultrasound in Medicine &Biology, 2010, 36(1): 111-129.

[19]孙 正, 杨 宇. 一种基于snake模型的IVUS图像序列三维分割方法[J]. 工程图学学报, 2011, 22(6):25-32.