冠脉造影图像是临床诊断冠心病的“金标准”。其成像过程中，受心脏跳动、成像设备等因素影响，图像不可避免地会引入一些噪声，使得图像中的血管脉络一定程度上模糊不清,尤其是血管边缘与细节的模糊，影响医生对病情的诊断。

近年来，不少学者致力于研究冠脉造影图像中血管脉络特征的提取。常见的提取冠脉造影图像中的血管信息方法主要有数学形态学方法^{[1, 2, 3]}、匹配滤波方法^{[4, 5, 6]}、基于中心线的方法^[7]等。其中，数学形态学方法是通过膨胀和腐蚀2种最主要的形态学操作实现图像中目标的扩张或收缩，以先腐蚀后膨胀的开操作和先膨胀后腐蚀的闭操作作为基础，实现图像中的边缘检测、图像增强和特征提取等相关操作；匹配滤波方法是将多个线性滤波器与卷积运算相结合的方式提取图像中的目标信息，通常还需要其他相关图像处理操作才能得到图像中的血管信息；基于中心线方法的基本思想是提取出图像中整个血管树的中心线，然后通过测量血管在某位置的半径重构图像中的血管脉络。但由于血管的形状复杂，血管的直径和造影剂的密度也会发生变化等，使得现有的方法并不能有效且完整地将冠脉造影图像中的血管信息提取出来。

经验模态分解是一种新的非平稳信号分析处理方法，可用于冠状动脉造影图像的图像处理过程^{[8, 9, 10]}。本文提出一种基于经验模态分解的血管特征提取方法：采用经验模态分解将造影图像分解成若干图像层，形成图像层组，该组中有的图像层包含了大量有效的血管信息而有的则包含的是肋骨、肺或其他组织等形成的非血管信息，丢弃含非血管信息的图像层，对含有血管信息的图像层进行血管特征提取，依据连通区域一系列属性提取出图像层中清晰明显的血管段，将各图像层提取出的血管信息相互参照并重构血管脉络明显的心血管图像。此外，本方法对冠脉造影图像经验模态分解之前须进行图像预处理操作，这是因为经验模态分解方法的基本原理是根据目标的变化频率而将图像中信息逐层分离出来，但受现有投影设备的限制，在冠脉造影图像的四周存在从有像素区到黑色无像素区的过渡区域，该区域的频率变化最大，如果不进行图像预处理进行图像像素填充使该区域实现平滑过渡，会在经验模态分解结果中出现严重边缘效应。

在对冠脉造影图像进行特征提取之前，须对图像中无像素黑色区域进行像素填充，使该区域像素点的灰度值实现平滑过渡，以抑制此后经验模态分解过程中边缘效应的出现。

受投影设备的限制，冠脉造影图像的四周会出现黑色区域，该区域不包含任何血管信息且频率变化最大，由于经验模态分解是根据目标的变化频率而将图像中信息逐层分离出来，如果不对此区域进行像素填充则会在分解结果中出现严重边缘效应。图 1所示的是临床常见的冠脉造影图像。观察图像可知在黑色区域与圆形亮色区域的相接处，像素点的灰度值变化最为剧烈，在图像经验模态分解过程中此处会出现伪血管信息。

图 1 冠脉造影图像Fig. 1 Coronary angiography images

图选项

为此，本文对图像中无像素的黑色区域进行像素填充，填充的过程如下：

1)将图 1的冠脉造影图像划分为如图 2所示的4个区域，箭头表示该区的像素填充方向。每个区域的填充方向各不相同，依次为左上、右上、左下和右下。

图 2 冠脉造影图像分区Fig. 2 Partition of angiography images

图选项

2)定义像素填充模板。本文是在待处理像素点的邻近位置选取大小合适的矩形填充模板，由此模板计算得出待填充点的像素值。例如，对填充方向为左上方的1区来说，采取的操作是在待填充点的右下方取一个大小为m×n的矩阵。根据式(1)计算该矩阵的平均灰度值，将该值作为待填充点的灰度值。如图 3所示，P是待填充点，本文中取m=4,n=4。

图 3 4×4大小的矩阵模块Fig. 3 Partition’s 4×4 matrix module

图选项

3)设定阈值并通过阈值T将待填充区域内的所有点分类成待填充点和保留点。如果一个像素点的灰度级小于阈值T则使用填充模板进行灰度级填充，如果大于阈值T则保留该像素点的原灰度级，不做灰度级填充处理。阈值T的设定方法是：对于选定的图像填充区域，首先统计该区域内每个灰度级像素的出现概率，即得到该区域的灰度直方图，然后概率小于某一数值的所有灰度级中的最大值作为阈值。本文实验中概率取1.16×10^-2，约为0.011 6。

4)依据图像填充模对图像①区由下至上、从右到左扫描，以填充黑色区域中的所有待填充点。对其他待填充区域进行类似处理，唯一不同的是矩阵模块与待填充点的邻接方向发生了变化。

5)图像二次填充以去除毛刺。经过4)对冠脉造影图像进行初次扫描及像素填充之后，部分像素点的过渡处仍存在小的毛刺，为了更好地实现图像像素平滑过渡的效果，重复4)对图像进行二次扫描填充。冠脉造影图像经像素填充预处理后的图像如图 4所示。下文中提到的冠脉造影图像均指预处理后的图像。

图 4 预处理后图像Fig. 4 The result of image preprocessing

图选项

本文采用经验模态分解方法将造影图像分解成若干图像层，利用了血管以下特性：在三维空间中为管状结构，经X光照射后形成冠脉造影图像中的血管段，其在血流方向上灰度级变化缓慢，而在垂直于血管血流方向上灰度级变化剧烈。如图 5所示，图像灰度级在λ₁方向上变化较平缓，在λ₂方向上变化剧烈。

图 5 心血管的空间结构Fig. 5 Spatial structure of the cardiovascular

图选项

以图 4为实验对象，图像通过经验模态分解后的结果是原图像被分解成不同的图像层，每一层中包含图像中不同变化频率血管特征信息，图像分解结果如图 6所示。通过对照原冠脉造影图像与图像分解后的图像层组发现，图像中的血管信息主要分布在imf₁、imf₂和imf₃中。因此，在对各分解层信号的后续处理过程中，本文只针对前3层进行处理。

图 6 图像EMD分解Fig. 6 Image EMD decomposition

图选项

针对冠脉造影图像强背景噪音的特点，经验模态分解虽然能很好地将造影图像中的血管信息和背景信息分离出来，并能将图像中的血管根据变化频率的不同进行分级，生成含有不用尺度级血管的图像层组，但是经验模态分解算法在分解图像产生若干个图像层的过程中某些情况下可能发生模态混叠现象。此现象会造成图像层中的血管特征信息失真，使得图像层中的信息不能真实反映冠脉造影图像中血管脉络信息，所以采用小波变换模极大值方法实现图像层间的信息参照以消除模态混叠造成的信息失真现象。

对冠脉造影图像分解后的图像层进行层间信息参照，主要是imf₁、imf₂和imf₃之间信息的参照。以图 6中图像为例，得到如图 7所示的极大值序列。

图 7 各图像层中的极大值Fig. 7 The maxima of each image layer

图选项

针对imf₃中的极大值点，如果在图像层imf₂相同图像位置上，该图像位置周围一定的区域范围内检测到极大值点，则imf₃中的此类点属于血管特征信息，将该类点标记并生成标记点集合Mark₁。同样的参照方法对imf₂图像层中的极大值进行标记，以imf₁图像层中的极大值点为参照信息，对Mark=Mark₁∪Mark₂中的极大值点，如果在图像层imf₁相同图像位置上，周围一定的区域范围内检测到极大值点，则imf₂中的此类点为标记点，生成的标记点结果标记点集合Mark₂中。由此可知，造影图像血管信息层间的信息参照结果是一个标记点集合imf₁，其中Mark=Mark₁∪Mark₂。在图像重构过程中标记点集合即为结果图像中血管树生成过程的参照信息，包含标记点的连通区域即为血管脉络树中的特征信息，否则该连通区域属于背景信息，做丢弃处理。采用类似小波变换的模极大值思想，对imf图像层进行层间极大值参照处理，结果得到一个标记点集合Mark{m_i(x,y)|i=1,2,…,k,1≤x≤M and 1≤y≤N}。依据标记点集合进行图像血管脉络重构。

依据小波变换模极大值的思想对各图像层间的血管信息进行层间参照，得到标记点集合。在对各层图像进行信息重构以实现血管脉络树重构时，参照标记点集合，实现各层信息的信息综合，从而得到最终的血管树，获取冠脉造影图像中的完整血管信息。算法1详细描述了本文中冠脉造影图像中的完整血管信息的提取过程。

算法1    基于EMD和模极大值的血管提取算法 输入：插值预处理后图像；

输出：含有血管脉络树的结果图像。

1)插值预处理后图像进行条件初始化: r₀(t)=x(t),i=1；

2)whiler_i(t)≥2 // r_i(t)是x(t)的“残余量”

3)do { //提取第i个本征模态函数；

4)    条件进一步初始化: h₀(t)=r_i(t),k=1；

5)        do{

6)            计算h_k－1(t)的局部极大值与局部极小值；

7)            利用三次样条对局部极大值与局部极小值序列进行插值；

8)            分别形成h_k－1(t)的上包络线和下包络线；

9)            计算h_k－1(t)的上、下包络线的均值，求得平均包络线m_k－1(t)；

10)            定义h_k(t)=h_k－1(t)－m_k－1(t)；

11)            k=k+1；}

12)          while h_k(t)不满足本征模态函数IMF的标准；

13)          imf_i(t)=h_k(t); //生成本征模态函数

14)            r_i(t)=r_i－1(t)－imf_i(t)；

15)          i=i+1;}

16)    分解完成，生成imf图像层组合残差图像；

17)图像层组分成包含血管信息的有效层和包含背景信息的无效层；

18)有效图像层进行极大值检测，层间信息相互参照，生成标记点集合；

19)依据Mark集合重构血管脉络树；

20)生成结果图像，血管信息提取过程完成。

依据标记点集合对冠脉造影图像进行重构，将各个有效图像层中的血管信息进行综合，生成结果图像。图 8 给出图像重构过程中不参照Mark集合和参照Mark集合的结果图像效果对比。在图 8(a)中，结果图像是直接将经验模态分解后的有效图像层中提取出的血管信息直接相加得到，图 8(b)中的结果图像则在图像重构过程中，参照标记点集合，将包含标记点的连通区域保留在结果图像中，不包含标记点的连通区域看做噪音处理。对比发现，通过模极大值进行图层间的信息相互参照，生成的标记点集合能很好地实现图像血管脉络树的重构并有效抑制噪音。

图 8 结果对比Fig. 8 Result contrast

图选项

为了更好地说明本文提出的算法能有效提取冠脉造影图像中的血管信息，将本文算法的血管提取结果与现有常用的Hessian矩阵、Gabor滤波器、直方图均衡化(HE)、线性反锐化掩模(LUM)等算法结果作比较。

图 9给出了Hessian矩阵、Gabor滤波器、直方图均衡化(HE)、线性反锐化掩模(LUM)和本文提出的基于EMD和模极大值的血管提取算法的实验结果对比。

图 9 算法对比Fig. 9 Algorithm effect contrast

图选项

以冠脉造影图像为例进行血管信息提取结果对比，图 9(b)给出了Hessian矩阵对造影图像的血管提取结果，该方法能将血管与背景对比度非常强的区域中的血管段提取出来，有很好的血管连续性，但是不足的是对于血管与背景对比度稍微较弱的区域则无能为力，能提取出的血管信息非常有限，提取效果不太理想。图 9(c)中给出了Gabor滤波器对造影图像中的血管信息进行过滤提取，能够将血管脉络的轮廓提取出来，但是得到的结果图像中血管特征并没有很好地凸显出来，血管树的连续性也较差。图 9(d)给出了造影图像直方图均衡化后得到的结果图像，该图像在一定程度上提高了冠脉造影图像的图像效果，但是对于背景噪音强、图像对比度低的图像区域还是存在背景噪音过强的缺点，整体来说还是能够将血管脉络结构较好地表达出来。图 9(e)给出了线性反锐化掩模对造影图像的处理结果，该方法能够实现造影图像中血管的细节和边缘突出表达，但是分段函数产生的量化误差较大时会将图像中的边缘区域对比度过度增强，产生伪影现象。图 9(f)给出了本文算法对造影图像的血管提取结果。

针对冠脉造影图像噪音大和血管边缘模糊等问题，提出的基于经验模态分解的血管特征提取方法，是将造影图像分解成若干图像层，形成图像层组，对图像层进行血管提取，将各层中有效的血管信息汇总在结果图像中，得到血管脉络完整的心血管图像。

对比其他算法的实验结果，文中算法得到的结果图像有极好的视觉效果，去除了造影图像中强背景噪音，使得图像更便于血管信息的直观分析和观察，同时提取出的血管树结构有很好的连续性和完整性，造影图像中血管特征信息，如血管方向、狭窄程序和血管分叉点都有很高的精确度，算法具有很好的实际医用价值。