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

张也，张勇德，沙宪政

中国医科大学 公共基础学院，沈阳市，110122

0 引言

视网膜是眼球壁内层的透明薄膜，血管从视盘向外延伸，由粗变细，血管灰度由大变小。SIMON等[1]提出视网膜血管形态具有高度唯一性和防伪性，可作为新的生物识别特征。TOWER教授[2]则在研究大量样本后得出结论，视网膜血管形态结构因人而异，即使是孪生子，视网膜血管结构也不相重。而且视网膜位于眼球最内层，与外界无接触，常态下形态结构不受外界因素影响，所以视网膜具有高度稳定性。同时视网膜血管是全身微循环系统的重要组成部分，在疾病预测、诊断方面也有很多应用。例如视网膜血管弯曲度改变是青光眼、高度近视、阻塞性睡眠呼吸暂停等疾病的表现[3]。视网膜血管检查与蛋白尿的检查及血压监测结合起来，可成为临床筛查慢性肾脏病的更有效方法[4]。因此，视网膜血管的准确分割，对疾病的检测和身份识别都具有重要价值。

目前，视网膜分割方法可分为传统方法与机器学习方法。传统方法有CHAUDHURI等[5]提出二维高斯匹配滤波器，该算法对大血管增强效果较好，分割速度快，但是对小血管增强效果较差。HOOVER等[6]利用血管局部信息，提出了多阈值探针血管检查方法去归类检测到的像素，算法假阳性率高，但是计算复杂、运行时间长。MENDONCA等[7]提出自适应阈值分算方法，分块并计算各部分的最佳阈值，再合并获得最终结果。陶治江等[8]改进了传统基于相位一致性的血管分割算法，初步增强图像后使用高帽变换进一步增强，最后使用形态学方法消除二值图像中的块状噪声。GUO在文献[9]中采用水平集方法，使水平集曲线逼近目标血管边缘。由于图像中噪声的存在，该方法分割精度不高。

机器学习方法有MARIN等[10]应用了基于矩不变特征的神经网络；LUPASCU等[11]提出了一种基于特征的AdaBoost分类器的类来自动分割视网膜血管；ASLANI等[12]提出了一种基于混合特征集的视网膜分割算法，算法综合了其他不同算法的特征，构成了一个17维的混合特征向量。这些特征包括Gabor 滤波器、形态学高帽变换、COSFIRE 滤波器等。算法使用一个快速而简单的随机森林分类器来训练和分割血管。LUO[13]将 CNN 与条件随机场（CRF）串联，提出了CNN-CRF 网络模型用以分割血管。MANINIS[14]基于 CNN 提出了视网膜图像深度理解网络（DRIU），该网络可以同时分割视盘和血管像素。

传统算法受图像噪声影响较大，细小血管分割效果差，细节丢失较为严重。机器学习方法需要大量的训练样本，运算速度慢，算法向硬件移植较为困难。因此本研究针对上述两类算法的缺陷，对传统的高斯匹配滤波算法进行了优化，设计了多尺度高斯匹配滤波器，并将不同尺度的滤波结果进行叠加，然后采用二维最大信息熵阈值分割算法对增强后的图像做二值化处理。

1 视网膜血管分割

1.1 视网膜图像预处理

本算法在DRIVE数据库中进行试验，该数据库中图像全为彩色图像，在进行后续的图像分割、图像识别、图像分析等上层操作之前，需要先将其转化为灰度图像。我们采用通道分离方法对彩色图像进行灰度化处理，结果如图1所示。与红色通道图像图1（a）与蓝色通道图1（b）的灰度图相比，绿色通道图像图1（c）亮度适中、灰度值分布均匀、血管与背景对比度强。因此选择绿色通道图像进行分割。

图1 视网膜图像预处理Fig.1 Preprocessing of retinal image

1.2 图像增强

视网膜血管横截面灰度强度呈倒高斯曲线分布，血管中心线处灰度值低，边缘灰度值高。CHAUDHURI等[5]设计了高斯匹配滤波器，用其分布来模拟血管横截面灰度强度分布，对血管分段进行滤波，当一段血管长度小于滤波器长度时，可将该血管段近似看做直线。若血管宽度与高斯核的尺度相匹配，则滤波器输出值最大。

定义二维高斯核函数为：

其中σ是高斯核的宽度，l为高斯核的长度，血管从视盘中心出发，向多个方向延展，旋转高斯核以对多方向血管进行滤波。假设p(x,y)为核函数中一个离散点，旋转矩阵为：

θi(0≤θi≤π)为第i个核函数的角度，p(x,y)旋转后的坐标值为则第i个模板核函数为：

其中Z为模板领域，取值范围为：Z={(u,v),|u|≤3σ,|v|≤l/2}。

Chaudhuri算法的滤波核参数设置单一，分别为l=9，σ=2，12个方向。该算法简单、运行速度快，参数设置能有效分割大血管，但对小血管的分割效果差。该研究针对高斯匹配滤波响应对血管管径依赖的问题，对其进行了改进和优化。

使用大尺度高斯匹配滤波核对图像增强结果如图2所示，其中参数设置为l=10.8，σ=1.9，8个方向[15]。从图像中可以看出，该算法对粗血管分割效果较好，抗噪性强，但是对细小血管的分割效果差，且较细血管无法与背景区分或者容易造成血管的断裂。

图2 大尺度高斯匹配滤波器增强图像Fig.2 Enhanced image with large-scale Gaussian matched filter

使用小尺度的高斯匹配滤波核对图像增强结果如图3所示。经多次试验，其参数设置为：l=5，σ=0.1，18个方向。使用小尺度滤波核，可以有效增强图像中的细小血管。同时也增强了图像中的一些条状噪声，且对存在中央反射的粗大血管增强效果差。

图3 小尺度高斯匹配滤波器增强图像Fig.3 Enhanced image with small-scale Gaussian matched filter

如图4所示，（a）图中血管存在中央反射，中心线区域的灰度值较大，边缘部分的灰度值较小。使用小尺度滤波器对其增强后，如图（b）所示，增强后的大血管近似可看做两条并列细小血管，血管中央反射区的灰度值与边缘处的灰度值差距较大，会造成错分。

图4 小尺度滤波器对存在中央反射的血管的增强Fig.4 Enhancement of vessels with central reflection by small-scale filter

针对单尺度匹配滤波各自的缺陷，本算法使用了多尺度高斯匹配滤波器，将两个尺度滤波器相结合，定义大尺度高斯匹配滤波核的滤波响应为φ1，参数设置为l=10.8，σ=1.9，8个方向；定义小尺度高斯匹配滤波核的的滤波响应为φ2，参数设置为l=5，σ=0.1，18个方向；则总的滤波响应为：

滤波增强结果如图5所示。从图中可以看出，粗血管与细小血管都得到了有效增强，解决了传统高斯匹配算法对存在中央反射区血管的误分问题，前景目标的灰度值与背景灰度值对比度得到了增强。

图5 多尺度高斯匹配滤波器增强图像Fig.5 Enhanced image with multiscale Gaussian matched filter

1.3 图像分割

视网膜灰度图经多尺度滤波器增强后，目标与背景分别占据不同灰度级，因此本研究采取阈值分割来对图像做二值化处理。阈值分割法可分为基于一维直方图的阈值分割法和基于二维直方图的阈值分割法。常见的基于一维直方图的阈值分割法有最大类间方差法、最小交叉熵法、最大熵法等，但是这些方法没有利用到图像的局部空间信息，当目标和背景的分布相互重叠和图像中存在噪声干扰时，分割效果不佳，如图6所示，使用大津阈值分割法[16]分割图像后，血管存在断裂、误分等问题。

图6 大津阈值分割法分割结果图Fig.6 Image using OTSU threshold segmentation

本研究使用基于二维图像熵的阈值分割算法。图像熵是一种图像特征的统计形式，反映了图像中平均信息量的多少。图像最大熵阈值分割的原理为使选择的阈值分割图像目标区域、背景区域两部分灰度统计的信息量为最大。

图像二维直方图定义如下：

其中，M×N表示图像大小，ni, j表示图像中灰度值为i，邻域灰度平均值为j的像素个数。

二维直方图的平面示意图如图7所示：

图7 二维直方图的平面示意图Fig.7 Plane sketch of two-dimensional histogram

横轴i表示点灰度值，纵轴j表示邻域灰度值。其中区域C和区域B表示目标和背景像素，区域A和区域D表示边界和噪声信息。阈值向量(t,s)，t表示灰度值，s表示像素邻域均值（通常是8邻域）。

对于L个灰度级的图像，设在阈值(t,s)定义区域1和2的概率P1，P2：

定义二维离散熵H一般表示为：

对各区域概率进行归一化处理可得区域1的二维熵：

同理区域2的二维熵：

其中，H1、H2为：

那么整个图像中目标和背景熵之和的函数：

根据最大熵原则，存在最佳的阈值向量满足条件：

后处理：分割后的二值图像还存在一些噪声，采用形态学方法将面积低于100的连通域去除。

2 实验结果及分析

为测试该算法的有效性，本研究在国际上公开的DRIVE[17]数据库中进行了实验。DRIVE数据库中的数据分为测试集和训练集，每个数据集中包含20幅彩色眼底图像和与之对应的专业医师手工分割的视网膜血管图像。该图像库中的视网膜图像采用用Canon CR5非散瞳3CCD照相机拍摄，视场角为45o。图像分辨率为768×584，R、G、B三个分量各有8 bit，图像全部视觉区域（Overall Visual Field,FVO）的直径大约是540个像素点。

这两个算法是传统算法中具有代表性的算法，所以用它们与本算法进行比较。采用图像来对比分析本算法性能。图8是各算法对DRIVE数据库中第一幅图像的分割结果。图(a)为CHAUDHURI等[5]提出的基于高斯匹配滤波器的分割算法，图（b）为FRANGI等[18]提出的基于Hessian矩阵的视网膜分割算法，图(c)为专业医师手工分割图像，图（d）为本算法的分割结果。对比可看出，本算法分割结果细节更丰富。

图8 本算法与其他分割算法的比较Fig.8 Comparisons between the proposed algorithm and other segmentation algorithms

为进一步量化分析算法性能，本研究采用了灵敏度（sensitivity，Sn）、特异性（specificity，Sp）、准确度（accuracy，Acc）这三个评价指标。灵敏度代表正确分割血管的比例，灵敏度越高，表示对血管的正确识别率越高。特异性代表正确识别背景的比例，特异性越高，表示对背景的识别率越高。准确度表示正确识别血管和背景的比例。准确度值越大，代表正确识别率越高。

其计算公式如下：

TP（True Positive）为真阳性，表示正确分割的血管像素个数；TN（True Negative）为真阴性，表示正确分割的背景像素的个数；FP（False Positive）为假阳性，代表将背景像素分割为血管像素的个数；FN（False Negative）为假阴性，代表将血管像素分割为背景像素的个数。

本算法与其他算法的性能比较结果如表1所示，从表中可以看出，本算法的灵敏度较高，达到了0.803，准确度高于其他算法，而特异性略低于其他算法。

表1 不同算法的性能比较Tab.1 Performance comparison of different algorithms

3 结束语

视网膜血管的精准分割对身份识别、疾病诊断都具有重要意义。本研究针对单尺度的高斯匹配滤波算法的缺陷，设计了多尺度多方向高斯匹配滤波算法。实验结果表明，本算法可以分割出更多的细小血管，且在灵敏度这一性能指标上优于多数经典算法。本算法模型参数调优简单，分割速度快，具有较高的临床价值和身份识别应用价值。