基于融合相位特征的视网膜血管分割算法

李媛媛，蔡轶珩，高旭蓉

(北京工业大学 信息学部，北京 100124)(*通信作者电子邮箱caiyiheng@bjut.com)

0 引言

眼底图像具有获取简便、成像清晰等优势，可以帮助人们尽早地检测眼部疾病，及时治疗。视网膜血管作为眼底图像的重要组成部分，是唯一可以通过非入侵的方式观察到的微小血管，而且一些全身性的疾病，如糖尿病、高血压等，往往会引起视网膜血管发生病变，所以通过对眼底视网膜血管图像的分析来进行一些疾病的早期检测和辅助治疗是一种非常好的医疗手段。

鉴于视网膜血管分割的重要性，国内外学者提出了许多有效的视网膜血管分割算法，主要分成非监督(例如，血管追踪[1]、匹配滤波[2]、基于模型[3]等)和监督[4-9]两类。其中，监督的方法是将眼科专家手工标记的视网膜血管图像作为标准样本，并基于给定的血管特征通过训练学习优化分类器参数得到最终的分类器，然后将此分类器用于眼底视网膜血管的分割。由于监督的方法依据已有的先验标记信息来设计，所以基于监督方法的血管分割效果普遍优于基于非监督方法的血管分割效果；但是基于监督的分割方法需要根据提取的图像特征及标准标记信息构造大量训练样本，在训练学习的过程中逐步优化分类器参数，所以其训练时间较长，分割效率不高，因此，在基于监督的分割方法中，提取优质的图像特征用于训练和分割，不仅可以提高分割的准确率，而且可以显著减少训练样本数和训练时长，提高算法效率。

Soares等[4]提取眼底视网膜图像的Gabor变换和像素灰度特征，利用贝叶斯分类器实现血管检测。其中，不同尺度和不同角度的Gabor变换可以对不同宽度不同走向的血管进行检测，而且Gabor变换可以滤除噪声且对光照变化不敏感；但是，该方法只考虑了图像Gabor特征和灰度特征，对细小血管的检测不够准确。吴奎等[5]在Gabor特征的基础上加入组合线性算子提高了分割的准确率。吴奎[6]提取了眼底视网膜图像的组合位移滤波响应(Combination Of Shifted FIlter REsponses， COSFIRE)特征，采用贝叶斯模型来进行血管分割。其中，对称和非对称COSFIRE滤波器的使用可以检测出连续血管及血管末端，得到更优质的血管特征。slani等[7]提取眼底视网膜图像的Hessian矩阵、高帽变换、Gabor变换、条带选择组合位移滤波响应(Bar-selective Combination Of Shifted FIlter REsponses, B-COSFIRE)特征以及像素灰度特征用于分类。该方法中作者评估了每个特征在混合特征向量中的重要性，其中B-COSFIRE滤波和Hessian矩阵的预测重要性指标较高。Zhu等[8]提取眼底视网膜图像的形态学、相位一致性、Hessian矩阵、向量散场以及局部特征用于分类。其中，向量散场特征具有较高的特异性，形态学特征有较高的敏感性。朱承璋等[9]提取眼底视网图像的不变矩、灰度共生矩阵、高斯拉普拉斯(Laplace Of Gaussian， LOG)算子结合高斯二阶导、相位一致性及Hessian矩阵特征用于分类。其中，LOG算子结合高斯二阶导特征有较高的敏感性，Hessian特征有较高的特异性，相位一致性特征有较高的准确性；但由于相位一致性特征主要用于血管边缘检测，提取出的相位一致性特征中血管中心部分未能检测出来，所以其敏感性低。

上述基于监督的眼底视网膜血管分割算法大多提取较多的图像特征，对每个像素构造高维的特征向量用于分类器训练与分割，耗时长。此外，现有文献中大多直接使用眼底视网膜图像的相位一致性特征。相位一致性特征虽对血管边缘信息有很好的检测效果，但其对血管中心部位检测不足，不适用于基于像素分类的视网膜血管分割算法。通过对常用的眼底视网膜图像特征的分析和研究，本文基于预测重要性指标较高的B-COSFIRE特征、Hessian矩阵特征以及对方向性检测敏感的Gabor变换特征，针对相位一致性特征对血管中心部位检测不足的问题，提出了基于融合相位特征的视网膜图像血管分割算法，该算法克服相位一致性特征对血管中心检测的不足，且提取较少的特征，实现高效的视网膜血管的分割。

1 基于融合相位特征的视网膜血管分割算法

本文算法是由预处理、特征提取、分类器分类以及后处理4个模块组成，具体算法流程如图1所示。预处理后，提取视网膜图像的Hessian矩阵特征、Gabor特征、B-COSFIRE特征以及相位特征，对图像中每个像素构造一个4D的特征向量，采用支持向量机(Support Vector Machine, SVM)分类器进行训练与分类，实现了高效的视网膜血管分割。

图1 本文算法流程

1.1 预处理

本文通过预处理过程增强眼底视网膜图像中血管与背景的对比度，便于后续的血管提取。在原始的彩色视网膜图像中：如图2(a)，血管与背景部分对比度不高，且存在图像噪声，不适合进行后续处理。而各类研究表明，绿色通道分量图像的对比度高于其他两个通道分量，且完整地保存了视网膜图像中的血管结构，故本文选取绿色通道分量图像进行后续图像处理；图2(b)为提取的绿色通道分量的图像。为进一步增强视网膜图像的对比度，抑制噪声，本文采用对图像增强效果较好的受限制对比度自适应直方图均衡化(Contrast Limited Adaptive Histogram Equalization, CLAHE)[6]作图像预处理；图2(c)为对绿色通道分量进行CLAHE的结果。预处理后，眼底视网膜图像中血管与背景的对比度明显增强，更适合后续的处理。

图2 预处理过程结果

1.2 特征提取

1.2.1 Hessian矩阵

根据眼底视网膜图像中目标血管呈线性管状结构的特点，本文采用对线性结构敏感的Hessian矩阵来提取像特征。Hessian矩阵是一个由多元函数二阶偏导数构成的方阵，它的特征值和特征向量可以很好地描述视网膜血管类似树杈形状的拓扑结构，并且通过滤波函数将眼底视网膜图像中的噪声去除。Frangi等[10]最先将Hessian矩阵用于眼底视网膜血管增强，并提出线性滤波函数。

由于视网膜血管的直径存在变化，单一尺度的Hessian矩阵的血管检测效果不佳，故将Hessian矩阵的差分运算与高斯核函数结合，通过改变高斯函数的标准差，来获得不同尺度下的Hessian矩阵。根据高斯函数的卷积性质，空间尺度导数Fab可由输入图像I与高斯函数的二阶倒数卷积得到：

(1)

通过分析Hessian矩阵的各特征值与其对应的线性目标关系，Frangi等[10]提出在σ尺度下的线性滤波公式：

(2)

1.2.2 Gabor变换

在视网膜血管图像特征中，Gabor变换具有良好的方向和尺度选择特性，可以检测不同方向与不同尺寸的血管，而且Gabor变换对光照变化不敏感，可以克服眼底视网膜图像光照不均的缺点，因此本文选择Gabor变换进行图像特征的提取。Gabor变换即窗函数为高斯函数的加窗傅里叶变换，有时可归为小波变换的一种。本文中Gabor变换的定义同文献[4]，即分析小波为Gabor小波的连续小波变换。

2D Gabor变换为：

(3)

(4)

在眼底视网膜图像血管检测中，Gabor变换参数设置大多同文献[1]，ε=4，K0∈[0,3]，旋转角度θ以10°为步长在10°～170°内旋转，并取不同方向变换的最大值为尺度a下的输出响应。膨胀尺度a的选择应与血管相匹配，包括可能的血管宽度。经过实验，本文取尺度a∈[2,5]，并将各尺度下的Gabor变换结果进行叠加作为最后输出。

1.2.3 B-COSFIRE

B-COSFIRE滤波模型由Azzopardi等[11]提出，这种模型模拟复杂感受野细胞的工作机制，对血管类的带状结构有很好的检测效果。在眼底视网膜血管检测中，可以通过对称和非对称的B-COSFIRE滤波器来检测出连续血管和血管末端，提取更完整的视网膜血管特征。B-COSFIRE滤波响应就是由一组高斯差分(Different of Gaussian, DoG)滤波响应的乘积加权几何平均得到。

滤波公式为：

(5)

B-COSFIRE滤波器通过自动配置过程来实现其对管状结构的选择性，图3为B-COSFIRE滤波器自动配置模型。在滤波器自动配置过程中，给定一个标准差为σ的DoG滤波器，给定中心点，如图3(a)、(b)中点1，该滤波器会沿着同心圆进行滤波响应，其中响应最大的位置就是表征强度变化大的点，即感兴趣的点，如图3(a)中的点2、3、4、5，图3(b)中的点2、3。这些感兴趣点的数目取决于所考虑的同心圆的数目和指定的原型模型，并用(σi,ρi,φi)来表示这些感兴趣的点，S={σi,ρi,φi|i=1,2,…,n}为这些点的集合，其中：σi为最大的DoG滤波器标准差，(ρi,φi)为这些点相对于中心点的极坐标，n为DoG滤波器个数。

图3 B-COSFIRE滤波器配置模型

为了提高各个点位置的容错性，对DoG滤波响应进行模糊和移位操作。模糊操作即计算每一个位置极坐标邻域的DoG滤波器的最大加权阈值响应，权重为DoG滤波响应与高斯函数Gσ′(x′,y′)的系数乘积，σ′=σ0′+αρi，σ0、α为常数，ρi为与中心滤波器之间的距离。移位操作即将模糊后的DoG滤波响应向φi的相反方向移动距离ρi，模糊移位后的DoG滤波响应为：

yi′)Gσ′(x′,y′)}

(6)

其中：Δxi=-ρicosφi；Δyi=-ρisinφi；-3σ′≤x′，y′≤3σ′。

B-COSFIRE滤波响应公式为：

(7)

1.2.4 相位特征

1)相位一致性特征。

对于一个由方波或三角波组成的傅里叶级数，在方波的阶跃或三角波的波峰波谷处，其各次谐波分量的相位相同。故，相位一致性认为图像信号相位一致性最大处就是图像特征[12]。在眼底视网图像血管检测中，相位一致性特征不受光照强度和图像对比度的影响，且对于血管的边缘有很好的检测效果。

Kovesi[13]提出了二维相位一致性特征的计算公式：

(8)

其中：o表示方向的索引；n为尺度信息；To为噪声补偿；Ano为加权函数；ε为常数，避免分母为零。具体推导过程及参数设置见文献[13]。

2)融合相位特征。

如图4(a)所示，相位一致性特征对血管边缘包括细微的分支血管都有很好的检测效果，但是无法检测出血管中心部分。故单一的相位一致性特征不适用于基于像素分类的眼底视网膜血管分割算法。为此，本文提出了一种新的融合相位特征，即将分别提取的相位一致性特征与Hessian矩阵特征进行小波融合，融合后的相位特征不仅保留了良好的边缘信息，其血管中心部分也被填充为白色，更适用于基于像素分类的视网膜血管的检测。

在小波融合过程中，源图像被分解成高频部分和低频部分，其中，高频部分包括源图像的细节和纹理信息，低频部分包括源图像的概貌和平均信息，因此，小波融合更容易提取图像的结构和细节信息。本文小波融合的具体步骤如下：

步骤1 对两幅源图像，如图4(a)相位一致性特征、图4(b)Hessian矩阵特征，分别进行二层小波分解，得到各自的高频子图像系数与低频子图像系数。

步骤2 将两幅源图像的高频系数和低频系数按照加权平均的融合规则进行融合，得到融合后的高频系数与低频系数。

步骤3 对融合后得到的系数进行小波重构，得到融合后图像，如图4(c)融合相位特征。融合相位特征即保留了相位一致性特征的良好血管边缘信息，且血管中间部分被填充为白色。

图4 相位特征融合过程

本文中提取的4D特征如图5所示。对视网膜图像中的每一个像素都构造了4维的特征向量，用于后续的分类。

1.3 SVM分类

本文使用公共的用于血管提取的数字视网膜图像(Digital Retinal Images for Vessel Extraction, DRIVE)数据库中的训练与测试集合，通过LIBSVM工具箱进行SVM训练和测试。

图5 本文方法的各特征图

训练过程 对DRIVE数据库的训练集中的每幅视网膜图像的每个像素构造4维的特征向量，并分别提取1 500个血管点和非血管点，归一化后，作为训练样本。通过LIBSVM工具箱对SVM进行训练，得到SVM模型。训练过程中选择径向基函数(Radial Basis Function, RBF)核函数作为核函数。为节省训练时间，本文中惩罚系数C选择常用默认值1，gamma选择0.01，训练得到SVM模型。

测试过程 从DRIVE数据库的测试集中任意选择视网膜图像并对其像素点提取特征向量，归一化后，作为测试样本，输入训练后得到的SVM模型，输出二值化的视网膜血管结构图。

1.4 后处理

采用连通域度量的方法，将输出图像中连通域面积小于25像素的孤立噪声点去除。

2 实验结果与分析

2.1 评价指标

为了对算法的分割结果进行评价，将本文的实验结果与DRIVE数据库中专家手工分割的结果进行比较，使用准确率(Accuracy, Acc)和受试者特性工作曲线面积(Area Under ROC Curve, AUC)对算法性能进行定量分析。

准确率是评价对视网膜图像中血管点和非血管点的正确分割的概率，计算公式为：

(9)

其中：TP(True Positive)为真正类，TN(True Negative)为真负类，FP(False Positive)为假正类，FN(False Negative)为假负类。

受试者工作曲线(Receiver Operating Characteristic, ROC)是反映敏感性和特异性连续变量的综合指标，ROC曲线以真阳性率为纵坐标，假阳性率为横坐标。受试者特性工作曲线面积(AUC)可由梯形规则计算得到，其值越接近1，表明算法的性能就越好。

2.2 实验结果

本文以DRIVE数据库中第1位专家手工分割结果为评价标准，实验结果如表1所示，本文算法进行血管分割所得平均Acc为0.957 4，平均AUC为0.970 2，结果较好。图6为本文方法在DRIVE数据库所有图像中分割准确率最高的最好分割结果，图中大量血管均被分割出来，可用于临床参考。图7为本文方法在DRIVE数据库所有图像中分割准确率最低的最坏分割结果，图中血管主干部分被分割出来，在发生病变的区域，病变部分被错误检测成血管。虽然图7中出现成块的错误血管区域，但这依然可以提醒医生此处异常，为医生进行疾病检测节省时间。

表1 本文算法分割结果

图6 本文方法在DRIVE数据库中分割准确率最高的最好分割结果

图7 本文方法在DRIVE数据库中分割准确率最低的最坏分割结果

为了分析融合相位特征用于像素分类的血管检测的效果，本文进行了不同特征组合的实验，结果如表2所示。组合1提取了眼底视网膜图像的Hessian矩阵、B-COSFIRE滤波和Gabor变换特征，构造3维特征向量用于分类器分割，所得Acc和AUC分别为0.957 3、0.965 8，该特征组合分割效果好。组合2在上述3维特征组合中加入相位一致性特征，构造4维的特征向量用于分类，但血管分割的Acc反而降低。因为相位一致性特征虽对血管边缘检测效果好，但其对血管中心部位检测不足，直接使用相位一致性特征反而降低了分割的准确率。组合3则在上述3维特征组合中加入融合相位特征，即本文方法，实验得到的血管分割的Acc和AUC均有所提升。结果表明，在组合特征中，融合相位特征比相位一致性特征更适用于基于像素的视网膜血管分割。

表2 不同特征组合分割结果比较

此外，为进一步验证融合相位特征的有效性，本文将提取的单一的融合相位特征与相位一致性特征，分别通过SVM进行分类，结果如表3所示，单一的相位一致性特征用于像素分类提取出血管的Acc和AUC分别为0.919 1、0.935 9，单一融合相位特征用于像素分类提取出血管的Acc和AUC分别为0.947 8和0.957 8，均高于相位一致性分割结果。实验结果表明，融合的相位特征比相位一致性特征更适用于基于像素分割的视网膜血管分割算法。

表3 单一特征分割结果比较

为了分析本文算法的有效性，将本文算法与近期基于特征提取且算法准确率较优的文献算法[6-9,11]，以及近期较热且具有相似复杂度的基于神经网络的文献算法[14-15]进行比较。表4为本文算法与文献[6-9,11,14-15]算法在DRIVE数据库上分割准确度的比较，图8为本文算法与文献[6,9]算法在DRIVE数据库上分割结果图的比较，图9为本文算法与文献[7,14]算法在DRIVE数据库上分割结果图的比较。如表4所示，与文献[6-9,11,14-15]算法相比，本文算法分割所得到的平均Acc与AUC总体较高。文献[9]算法的Acc较高，但其AUC低于本文算法，且如图8所示，相较于文献[9]分割得到的血管图像，本文算法分割得到更多的细小血管且血管连续性较好。文献[14]算法的AUC较高，但其Acc略低，且如图9所示，相较于文献[14]，本文分割虽然部分血管末端连续性不够，但本文算法分割出更多的细小血管。故从总体指标来看，本文算法结果较优。此外，相较于文献[8]算法构造的39维特征向量、文献[9]算法构造的23维特征向量，文献[14]算法采用的完全卷积神经网络，本文算法仅提取了4维特征通过SVM分类，降低了算法的复杂度，提高了分割效率。

表4 不同算法分割结果比较

3 结语

针对相位一致性特征对血管中心检测不足的问题，本文提出基于融合相位特征的视网膜血管分割算法。不同于文献[8-9]算法，本文算法没有直接使用相位一致性特征，而是将相位一致性特征同Hessian矩阵特征进行小波融合，得到一种新的融合相位特征。该特征既保留了相位一致性特征中良好的血管边缘信息，又克服了相位一致性特征中血管中心部位与背景相似的缺点。由于基于像素分类的视网膜血管分割算法需要对每一个像素的特征进行识别分类，因此，本文中的融合相位特征比相位一致性特征更适用于基于像素分类的血管分割算法。此外，本文对每一个像素提取Hessian矩阵、Gabor小波、B-COSFIRE、相位特征，构造4维特征向量用于SVM分类，提取较少的特征与样本，实现快速、高效的眼底视网膜血管的分割；但是，与现有大多的眼底视网膜血管分割算法一样，本文对于眼底视网膜中毛细血管的分割不足，损失了一些细小血管的信息。在接下来的工作中，如何对细微小的血管进行更高准确度的分割仍需要继续研究。

图8 文献[6,9]算法与本文算法结果比较

图9 文献[7,14]算法与本文算法结果比较