基于多尺度滤波的视网膜血管分割算法

梁礼明 盛校棋 蓝智敏 钱艳群

(江西理工大学电气工程与自动化学院 江西 赣州 341000)

0 引 言

血管是视网膜最重要的组成部分之一，视网膜血管分割和血管形态属性的划分，如长度、宽度、迂曲和角度，可用于各种心血管和眼科疾病的诊断、筛选、治疗和评估[1]，如糖尿病、青光眼[2]、硬性渗出物[3]和高血压[4]等全身性疾病。然而，通过眼科专家分割视网膜血管存在时效性差，且缺乏准确性。因此，研究血管自动分割和自动识别血管形态属性的算法成为辅助医疗诊断至关重要的一步。

现有的视网膜分割方法大致可划分为无监督与有监督方法。其中，无监督学习视网膜血管分割方法可以分为匹配滤波、数学形态学、血管追踪以及聚类等。文献[5]提出一种改进的Top-Hat变换应用于视网膜血管分割，使用不同半径的圆形结构元素来检测不同宽度的血管，进一步提高微小血管的分割率，但对噪声较为敏感。文献[6]提出一种基于贝叶斯概率的血管追踪方法，结合血管的连通性与灰度级等信息拟合血管结构，该方法较好地解决血管交叉处分割断裂的问题，但对微小血管分割精度较低。文献[7]利用主成分分析(PCA)提取血管的特征，并利用阈值进行分割，存在微小血管分割低、血管交叉处分割断裂等现象。

有监督学习需要专家提供具有标记的金标准视网膜图像，使用基于局部或全局图像的一组特征来训练分类器，充当先验知识并指导训练。文献[8]提出了一种融合多特征的随机森林分类器，很好地解决了血管交叉处分割断裂的情况，但存在将视盘误分割为血管的现象。文献[9]融合相位特征的血管分割算法，较好地解决了血管相位一致性特征检测不足的问题，但仍存在微血管分割不足的问题。文献[10]将血管分割看作为二分类问题，提出一种基于支持向量机(SVM)的混合5D特征区分血管像素与非血管像素，对噪声有较强的适应性，解决了视盘被误分割为血管的问题，但存在微小血管易断裂问题。

因此，针对现有的方法的不足与挑战，本文提出一种基于多尺度滤波器与AdaBoost分类器相结合的有监督的视网膜分割算法。首先利用二维K-L变换综合分析图像的三通道颜色特征信息，解决单一通道存在伪影和噪音大的问题，并采用受限对比度直方图均衡化(Contrast Limited Adaptive Histogram Equalization, CLAHE)增强视网膜图像并且移除低频的噪音，利用双边滤波Retinex进行伪影矫正。其次，利用多尺度形态学滤波、多尺度匹配滤波、2D-Gabor滤波及Hessian的各项异性提取较宽和微小血管特征向量，将以上先验信息与血管特征作为训练样本并归一化。使用专家分割结果作为标签，通过AdaBoost进行分类，去掉非血管像素，得到初步分割结果。利用AdaBoost分类器与多尺度多类型滤波器相结合，有利于改善分割效果，每个滤波器以不同的方式响应图像中的不同像素。最后，利用连通域面积等操作去除病灶与伪影，获得最终的视网膜血管分割结果。总流程如图1所示。

图1 本文算法总流程图

1 数据集与图像预处理

1.1 数据集

本文采用STARE[11]和DRIVE[12]眼底图像数据集作为实验数据。其中：DRIVE数据集各有20张训练集和测试集；STARE数据集仅含20张眼底图像，但并没有划分训练集与测试集。故本文利用留一法，选择一幅图像作为测试集，并使用数据集中的其他图像训练分类器，通过改变测试图像重复该过程，直到数据集中所有图像被用于测试一次。

1.2 二维K-L变换

采用二维K-L变换将视网膜图像三频带信息转换到一个主分量空间，进而降低算法复杂度，并有效地解决了人工通过加权的方式得到灰度图像的问题，其数据变换矩阵定义如下：

(1)

式中：αk=(αkR,αkG,αkB)T是图像p的协方差矩阵的特征向量矩阵；k=(1,2,3)；p=(pR,pG,pB)T表示RGB三通道图像。其主分量信息需要将数据进行协方差对角化，故定义协方差矩阵C(i,j)为：

(2)

(a) 绿色通道 (b) 第一主分量 (c) CLAHE图2 K-L变换图像与对比度增强图像

1.3 伪影矫正

考虑到视网膜图像包含光照不均匀与伪影等现象，因此利用改进的双边滤波Retinex[13]模型。经CLAHE增强的视网膜图像Ic，可视为反射率图像Rz与照明图像L的两个分量分块乘法：

Ic=Rz×L

(3)

本文定义x作为输入图像Ic的像素值，反射率图像Rz(x)的像素值x可通过计算输入图像Ic(x)和照明图像L(x)的对数之差来获得，Rz(x)定义如下：

Rz(x)=log(Ic(x)+1)-log(L(x)+1)

(4)

L(x)定义为：

(5)

归一化因子M(x)定义为：

(6)

式中:参数φ表示用于测量输入图像像素Ic像素与相邻像素之间空间关系的窗口大小，设置窗口大小为3×3像素[13]。

空间关系g(,x)定义为：

(7)

相似性关系s(,x)定义为：

(8)

式中：d(·)为欧氏距离；参数σd与σr分别表示视网膜图像空间扩展与强度扩展。函数g(,x)通过欧式距离测量Ic与空间像素关系；函数s(,x)代表Ic与强度相似性关系。

该模型的双边滤波法平滑视网膜图像，保持小邻域内像素x的强度值彼此近似，并且滤波强度基本一致，同时归一化项M确保所有临近强度的权重之和为1，以此充分保留血管的纹理特征，如图3所示。

(a) 传统Retinex (b) 本文Retinex图3 视网膜图像Retinex矫正

由图3可看出：传统Retinex矫正虽然整体亮度得到提高，并且抹掉大部分视网膜图像的伪影，但仍然存在大面积灰度信息不均匀的光照影响。而本文Retinex算法较好地矫正了视网膜图像的伪影现象，同时保留血管纹理特征等结构信息，较好地解决后期血管分割将黄斑和视盘误分割为血管的现象。

2 特征提取

2.1 形态学滤波

利用多尺度形态学Top-Hot变换[14]，在图像整体增强的同时，提取视网膜血管的微血管特征信息。通过控制图像边缘梯度信息控制因子wi，调整相邻血管像素尺度的差值，降低视盘与黄斑等特征信息的干扰，提高视网膜图像微小血管的多尺度亮、暗细节特征。多尺度Top-Hot模型定义如下：

(9)

式中：k是视网膜图像细节增强因子；Ir是输入经双边滤波的Retinex增强的视网膜图像；fT是多尺度Top-Hot增强后图像;Dopi与Dcli分别为视网膜图像亮细节与暗细节特征。控制因子wi的值为：

(10)

式中:x、y为经预处理后的视网膜图像邻域像素值；gi max与gi min分别为gi的最大值与最小值；gi是视网膜图像膨胀与腐蚀之差。wi的变化主要由视网膜图像的梯度信息决定，经多尺度Top-Hot增强的图像，再利用不同尺度下形态学低帽[15]变换提取血管特征，共计6维特征向量。

2.2 匹配滤波

多尺度高斯匹配滤波器[16]可以用来模拟眼底视网膜血管的整体走向与灰度信息，解决单一尺度下血管宽度变化较大时，血管匹配不全的弊端，其定义如下：

(11)

尺度函数Pi,j(x,y)由两个不同的响应函数产生，定义如下：

Pi,j(x,y)=Ri(x,y)Rj(x,y)

(12)

多尺度高斯滤波器的响应函数Ri(x,y)定义如下：

Ri(x,y)=fi(x,y)×im(x,y)

(13)

式中：im(x,y)是归一化的Retinex图像，其定义如下:

(14)

在本文的多尺度高斯滤波器中，ts为常数项t和滤波器尺度s的乘积；P1,2(x,y)与P3,4(x,y)别用来提取微小血管和宽的血管，归一化时设定0.75

2.3 Frangi滤波

基于Hessian矩阵的Frangi滤波器[18]，具有可将视网膜宽度不同的血管突出显示，并且由二阶导数在两个不同尺度上加以分析。Hessian矩阵特征值之间的差异可用于进一步的血管对比度增强和抑制非血管结构。基于Frangi滤波器的Hessian矩阵F(x)定义如下：

F(x)=maxδf(x,δ)

(15)

式中：感兴趣像素区域由x定义；δ是用于计算图像的高斯导数标准差，即血管的尺度方向；f(·)表示滤波器。原函数与高斯二阶偏导的卷积可降低二阶导数对图像噪音的敏感度。Hessian矩阵定义如下:

(16)

式中：Ir为经Retinex后的图像；Gδ(x,y)为高斯函数；Hessian矩阵中的Hxx、Hxy、Hyx和Hyy表示带有高斯核函数的二阶方向导数。定义λ1与λ2为H(x,y)特征值，由0<|λ1|≤|λ2|确定是否为血管像素和明显的灰度值变化，同时增强血管内部亮度，矫正血管中心亮线。为此构建血管响应函数Vδ(x,y)以区分背景与血管：

(17)

2.4 2D-Gabor滤波

2D-Gabor[19]滤波使用由正弦平面波调制的高斯核函数，其Gabor核具有不同权重的平行带，其核参数控制条带的大小、方向和位置，能够捕捉到血管空间尺度、位置和相位等局部结构信息，增强血管与背景对比度，且对图像伪影信息具有良好的适应性。本文给出一种基于小波变换的2D-Gabor滤波，设计如下：

首先，利用带有变换小波ψb,θ,a的标量积来定义连续小波变换Tψ(b,θ,a)：

(18)

式中：f为视网膜图像平方可积函数；x为图像像素值；参数b、θ、a分别表示位移矢量、旋转角、尺度；r为相位偏移量；ψ*为ψ的复共轭；Cψ与ψ分别为归一化参数和小波分析。

然后，给出2D-Gabor小波的定义如下：

(19)

(20)

其中:θ在Gabor小波变换的跨度以10°为步长[19]，从0°变化到170°。取每个尺度滤波结果的最大值作为该点的血管特征，该特征能有效区分背景与血管之间的灰度信息，共计12维特征向量。

以上4个滤波器特征图如图4所示。由于各个滤波器提取的维数特征不一，故将各个特征向量进行归一化到[-1,1]，并对每幅DRIVE和STARE的训练集选取7 500个血管与背景像素点，作为正反样本点，并将血管像素记为1，背景像素记为0。

(a) 形态学滤波 (b) 匹配滤波 (c) Frangi滤波 (d) Gabor滤波图4 多尺度滤波特征图

3 AdaBoost分类器及后处理

3.1 AdaBoost算法步骤

AdaBoost的学习过程是基于定义一个分类器，该分类器最小化训练集X的预测误差，这个最小化可以视为在训练集上沿着一组权重Di=(h1,h2,…,hT)(t=1,2,…,T)，最小化分类错误，具体步骤如下：

Stept1将经多尺度滤波器提取的视网膜血管特征向量中血管像素与非血管像素以1∶4的比例因子，作为训练数据集X=(x1,y1),(x2,y2),…,(xm,ym),yi∈{-1,1}，i=1,2,…,m,m为训练集数目，设定迭代次数为Tmax。

Stept4更新训练数据集的权值分布，提升样本在弱分类器中训练集错误分类的权重。

Stept5经过多次迭代，当所有样本错误率小于50%或者几乎为0，分类器的精度将不再提升，弱分类器通过检查所有特征的阈值来搜索强分类器H(x)，最终的强分类器由多个弱分类器加权组合而成，并将阈值应用于分类器的输出以决定分割的血管与背景像素。H(x)定义如下：

(21)

权重迭代规则定义如下：

(22)

3.2 后处理

AdaBoost分类器在根据数据集特征分割血管时，存在将块状小斑点区域误识别为血管。因此，本文构建一个8连通域的几何形态算子。该算子的阈值由图像的灰度统计直方图确定，通过统计去除小于30个像素组成的非细长和非连通区域，得到最终分割结果。由连通域面积(Area)和连通域宽(W)、高(H)信息构建如下几何算子[7]：

(2)W×H<3.5Area；

(3)Area<30。

4 AdaBoost分类器及后处理

本实验仿真平台为MATLAB R2018a，计算机配置为Intel(R)CoreTMi5- 7300 CPU@2.5 GHz,16.00 GB内存,64-bit Windows 10。

4.1 本文方法血管分割图像比较

图5展示了本文算法在健康视网膜图像与病变视网膜图像的分割结果。其中图5(a)与(e)为病变视网膜图像与健康视网膜图像灰度图，(b)与(f)为第二专家手工分割结果。观察(c)与(g)可知本算法在病变视网膜上会分割出多余的病灶与伪影区域，而在健康视网膜的分割结果比较理想，血管结构清晰，基本没有多余的类病灶区域的出现。(d)与(h)为由连通域除去小于30像素之后图像。通过观察(d)利用后处理的方法，分割结果得到了显著的改善，从而证明了本文后处理的有效性与必要性。

(a) Image0005 (b) 金标准 (c) 初分割图 (d) 后处理图

(e) image0163 (f) 金标准 (g) 初分割图(h) 后处理图图5 本文算法分割结果

4.2 不同文献血管分割图像比较

图6展示的是本文算法、文献[8]和文献[20]在DRIVE与STARE数据库的部分图像的分割结果。其中：(a)为第二专家金标准图像；(b)为本文算法分割结果；(c)为文献[20]利用全连接神经网络与条件随机场的血管分割图像；(d)为文献[7]融合多特征的随机森林血管分割算法。

(a) 金标准 (b) 本文算法 (c) 文献[20] (d) 文献[8]图6 不同算法分割结果

由图6前两行DRIVE数据集和第三行STARE数据集分割结果可知，文献[8]和文献[20]均出现视盘误分割为血管的现象，同时文献[20]丢失了大部份微血管信息，文献[8]血管出现过粗的现象。然而本文算法较好地解决了上述两种算法存在的缺陷，在视盘处基本不存在误分割的可能，假阳性较低，从而说明引入伪影矫正和二维K-L变化很好抑制了部分噪音信息，提高算法总体鲁棒性。并且本文相比上述算法能够分割处更多的微血管信息，泛化性较强。观察第二行分割结果可以看出，本文分割结果较文献[8]和文献[20]噪音含量极少，血管纹理结构清晰。

综上所述，本文算法相比以上算法具有较小的误分割率，能很好地减少视盘被误分割可能，在血管交叉处不易断裂且相邻血管不易相连，对伪影等病灶处有较好的鲁棒性。对于主血管末端的微小血管，本文算法能得到较好的保留，且不易存在微小血管断裂与分割过粗的现象，有较好的连通性与完整性，在保持精度与特异性高的情况下，尽可能分割出较多的血管。

4.3 视网膜分割评价指标

为了进一步定量分析本文算法对视网膜血管提取的算法性能，定义如下几个评价指标：

(23)

(24)

(25)

式中：TP、TN、FP、FN分别表示真阳性、真阴性、假阳性、假阴性。敏感度(Sensitivity)又称真阳性率，表示正确分类血管像素占真实血管像素的百分比；特异性(Specificity)表示正确分类的非血管像素占真实非血管像素的百分比; 准确率(Accuracy)表示正确分类血管和非血管像素占整个图像总像素的百分率。

4.3.1融合多种滤波器与单类型滤波器的比较

为了验证本文将AdaBoost分类器与多尺度多类型滤波器结合的有效性，给出在DRIVE数据库下单一滤波器取得的精度与灵敏度指标，并与本文算法的最终结果进行比较，如表1所示。

表1 各个滤波器与多类型滤波器性能比较

由表1可知单滤波器在AdaBoost下也可以得到较好的性能，但单一滤波器存在弊端，将其融合可以减弱单一滤波器相关缺陷。首先，通过基Hessian矩阵的Frangi可增强血管轮廓，多尺度形态学滤波增强微血管信息，降低粗血管与微小血管的对比度信息；2D-Gabor滤波能捕捉血管空间尺度，获得其相位信息；多尺度高斯滤波可以整合血管整体走向和灰度信息。然后将初步提取的血管特征信息通过AdaBoost多个弱分类器进一步与金标准比对提取血管特征，并加权组合得到强分类器完成最后血管的分类。该过程将无监督学习与有监督学习进行有机融合，使得总体算法精度从平均94%左右达到96%，灵敏度从平均76%提升到83%，更进一步说明本文算法融合4种多尺度滤波器比单一滤波器在AdaBoost分类器分割结果的有效性。

4.3.2不同视网膜分割算法性能比较

表2和表3给出不同算法与本文视网膜分割算法的性能，其中加粗的数值为分割结果最优异的值。选取DRIVE与STARE数据集中准确率、灵敏度和特异性三个指标作为不同算法间的比较标准。

表2 DRIVE数据集视网膜血管分割结果

表3 STARE数据库视网膜血管分割结果

由表2可知，现有算法在DRIVE数据集血管分割准确率均低于本文算法，说明本文算法在改数据集中具有较强的分割处血管与背景区域的能力。文献[7]利用PCA的方法分割血管，其特异性高于本文0.003 1，但其灵敏度与准确率分别低于本文算法0.016 8和0.249，说明文献[7]虽然将阳性较低，但存在血管分割不全和微血管分割严重不足的现象。文献[8]利用随机森林的血管分割算法，其灵敏度与本文算法相近，但文献[8]特异性低于本文算法0.018 8，说明本文算法误分割率远低于文献[8]。

从表3可知，本文算法在准确率和特异性方面高于现有算法。文献[10]灵敏度仅高于本文算法0.010 3，但另外两个指标低于本文算法，同时文献[10]采用的SVM分割算法还面临核函数选择和易过拟合的问题。文献[8]融合多种像素特征并利用随机森林算法分割血管，其灵敏度与准确率虽与本文算法较相近，但均低于本文算法0.004 9和0.002 4，其中特异性低于本文算法0.011 5，从而说明本文算法能在保证较高准确率与灵敏度的情况下，分割的血管图像假阳性率更低。

综上所述，本文算法相比以上算法具有较小的误分割率，能很好地减少视盘被误分割可能，在血管交叉处不易断裂且相邻血管不易相连，对伪影等病灶处有较好的鲁棒性。对于主血管末端的微小血管，本文算法能得到较好的保留，且不易存在微小血管断裂与分割过粗的现象，有较好的连通性与完整性，在保持精度与特异性高的情况下，尽可能分割出较多的血管，提高其应用价值。

5 结 语

本文利用二维K-L变换提取图像单通道信息，解决选取通道信息经验化的问题；利用双边滤波Retinex算法矫正图像的伪影等问题，较好地解决了大部分现存算法对存在伪影严重的视网膜图像分割不足的问题，并且进一步抑制视盘区域的对比度，解决后期被误分割为血管像素的难题。通过巧妙地选取具有不同属性的多尺度滤波器与AdaBoost相结合，充分发挥二者的优势，较好地解决了微血管分割不足的难题，并且对于存在病灶眼底视网膜图像中也具有较高的鲁棒性。考虑到数据集之间的不同图像具有不同的特征，本文通过组合这些滤波器，有效地克服单一特征或者滤波器存在的微血管分割不足与血管连通性不佳等问题，还可以进一步使分割方法更加稳健。同时，本文算法仍存在部分微血管分割断裂的现象，故下一步主要工作能有效重构微血管信息增强算法拓扑结构。