基于Hesse矩阵和多尺度分析的视网膜动静脉血管管径测量方法

肖志涛 崔 宁 吴 骏 耿 磊 张 芳 温 佳 童 军 刘晓婷 杨 嵩



基于Hesse矩阵和多尺度分析的视网膜动静脉血管管径测量方法

肖志涛 崔 宁 吴 骏*耿 磊 张 芳 温 佳 童 军 刘晓婷 杨 嵩

(天津工业大学电子与信息工程学院 天津 300387)

许多全身性疾病会引起视网膜血管管径及动静脉血管管径比例(Arteriolar-to-Venular diameter Ratios, AVR)的变化，因此对视网膜血管管径进行准确的量化分析对病情诊断具有重要的意义。该文提出一种视网膜动静脉血管管径及AVR的自动测量方法。首先，在分割血管网络的基础上，依据Hesse矩阵检测线状结构的优势，结合多尺度分析准确定位血管方向并计算血管管径；然后利用广义回归神经网络(General Regression Neural Network, GRNN)分类器对动静脉血管骨架线上的点进行准确分类；最后计算感兴趣区域(Region Of Interest, ROI)内的AVR。对REVIEW和DRIVE数据库进行实验，验证了所提方法的有效性。

视网膜血管；Hesse矩阵；管径测量；广义回归神经网络；视网膜血管管径与动静脉管径比

1 引言

视网膜血管管径与动静脉管径比(Arteriolar- to-Venular diameter Ratios, AVR)的改变可预测中风、糖尿病、高血压等全身性疾病的发病率[1]。眼科医生手动标记动静脉血管和测量血管管径准确率高，但耗时较大，不适合批量处理图像，因此研究自动准确的测量视网膜血管管径与动静脉血管比例的方法有助于提高医生的诊断效率。

正确分类动静脉血管是得到动静脉血管管径和AVR的前提，目前自动识别动静脉血管的方法大体分为有监督和无监督的分类方法。有监督的方法是在有标准图像的基础上提取特征训练分类器，例如文献[2]提出了训练K近邻分类器的方法；无监督的方法是在无标准图像的基础上设定判别准则，例如文献[3]以视盘为中心，将视网膜图像分成4象限，在每个象限内部提取血管特征，并应用模糊聚类对血管进行分类。目前现有的方法主要针对视网膜主干血管的分类，缺少对细小血管的研究。

现有血管管径测量的方法大体分为3大类。第1类为边界交点法，该类方法的关键是要准确得到血管方向。例如文献[4]提出了一种利用血管骨架的方向及空间位置，并采用改进的定向局部对比度方法测量管径的方法。第2类是横截面轮廓法，横跨血管的像素灰度近似高斯曲线，采用高斯曲线拟合横跨血管的像素灰度求得血管管径。例如文献[5]在基于视网膜血管的方向信息和拟合横跨血管像素点灰度的基础上，针对血管壁是否存在反射光的情况，拟合局部血管2维高斯曲线测量血管管径。第3类是建立模型法，是根据血管特征建立适当的模型测量血管管径，例如文献[6]提出了一种利用决策树和多分辨率的海尔曼模型的管径测量方法，拟合3维血管曲面模型计算血管管径。

目前现有的血管管径测量方法中边界交点法中存在不能测量血管分叉和分叉处管径的问题；横截面轮廓法中由于并不是所有的血管横截面灰度都呈高斯曲线分布，因而高斯曲线并不能完全拟合，管径测量会存在误差；模型建立法中有些方法需要设定参数不能达到自适应测量的效果或者不能达到全自动化的要求。因而当下迫切需要研究出一种精确且自动化测量血管管径的方法。

本文在分析现有动静脉血管分类和管径测量方法的基础上，结合Hesse矩阵和多尺度分析提出了一种新的自动测量血管管径的方法。根据Hesse矩阵检测局部线状结构的性质，结合Hesse矩阵和多尺度分析对眼底图像进行增强，并分割出完整的血管网络，本文方法的血管分割结果不仅能够保证血管的连续，而且既能够分割出较粗的主血管，又能够分割出毛细血管。在血管分割的基础上本文方法可以准确测量血管交叉和分叉处的管径，同时能够自动得到动静脉血管管径，以及感兴趣区域的AVR。通过对REVIEW和DRIVE数据库进行测试，结果表明本文方法具有较高的有效性。

2 基于Hesse矩阵的多尺度增强滤波器

视网膜血管为局部线状结构且宽度不一，采用对线状结构敏感的Hesse矩阵检测血管。对于2维图像，Hesse矩阵[7]是一个二阶方阵，可表示为

经过滤波后，视网膜图像中的血管会得到增强。为了从滤波结果中提取出血管，根据Hesse矩阵特征值与血管结构的关系，本文引入血管置信度函数来描述血管强度。单一尺度的血管置信度函数定义为[8]

图1 血管上Hesse矩阵的特征值对应的特征向量图

因此，为了描述不同尺度的血管，本文依据尺度空间理论，采取多尺度增强方式增强血管。当视网膜图像上通过迭代尺度()时会得到不同尺度下的，对于血管，只有当尺度因子与血管管径最匹配时值输出最大，取每个点处的最大值作为当前点的血管置信度。

3 基于Hesse矩阵的多尺度血管分割

3.1预处理

本文选取视网膜图像的对比度较高的绿色通道图像进行处理。视网膜图像中除血管呈局部线性特征之外，在视盘边缘处及背景中也存在部分线状特征，为了消除此类伪血管的干扰，首先要对绿色通道图像进行预处理。

为了降低伪血管特征对血管分割的影响，首先对绿色通道图像进行形态学平滑处理。形态学平滑处理是为了消除亮区域和暗区域的噪声，结果如图2(a)所示，可以看出去除了血管；令图2(a)减去绿色通道图像后反相，减弱了背景噪声突出了血管，结果如图2(b)所示；之后对图2(b)采用对比度受限的自适应直方图均衡化(Contrast Limited Adaptive Histogram Equalization, CLAHE)对视网膜图像进行处理，CLAHE同时具有自适应直方图均衡化和对比度限幅的优点，能提高视网膜图像血管的对比度和清晰度，如图2(c)所示。

3.2血管分割

在对图2(b)进行CLAHE增强后，进行基于Hesse矩阵的多尺度增强，结果如图3(a)所示。视网膜图像血管宽度由一个像素到十几个像素变化，为了增强各个尺度下的血管，根据文献[8]和实验仿真，将初始参数设为,，将范围设为0.5~10, 0.5为迭代步长。从图3(a)中可以看出，直接对图2(c)进行Hesse矩阵多尺度增强，不仅血管得到增强，同时也增强了部分背景噪声且噪声呈点状分布。为此在Hesse矩阵多尺度之前，采用各向异性耦合扩散方程[9]对图2(c)滤波以滤除背景噪声。由于该方程在血管边缘处滤波程度较小，在血管内部滤波程度较大，且滤波主要在图像边缘的切线方向进行，因此既能有效滤除噪声又能很好地保持血管的边缘信息，结果如图3(b)所示。对图3(b)进一步进行基于Hesse矩阵的多尺度融合滤波器增强，结果如图3(c)所示，可以看出血管边缘更加清晰，点状噪声得到抑制，残留的背景噪声呈块状分布。

为进一步突出血管并消除块状噪声，将图3(a)和图3(c)进行像素级相乘，结果如图3(d)所示，从图中可以看出血管结构得到增强，背景噪声则被进一步消除。最后用迭代阈值法将血管二值化，并进行形态学后处理，最终血管分割结果如图3(e)所示。

图2 预处理结果

图3 血管增强和分割结果

4 血管管径测量

本文在血管分割结果基础上测量血管管径。管径测量包括3步，首先提取血管骨架，然后沿垂直于血管的方向扫描找到血管边界点，计算两边界点之间的距离获得血管的像素级宽度。数学形态学细化方法提取血管骨架，结果如图4所示。

图6为视网膜图像中局部部分血管管径测量示意图，可以看出本文管径测量方法可以很好地作用于血管交叉和分叉处。

图4 血管骨架图

图5 血管管径测量示意图

图6 视网膜血管管径测量示意图

5 动静脉血管分类

为了分别得到动静脉血管管径并计算AVR，还进一步需要对动静脉血管进行分类。首先分析视网膜图像中的动静脉血管，动脉呈橘红色，颜色浅，较细；静脉呈暗红色，可见血管反光带，颜色深，较粗。动脉和静脉正常管径比例约为2:3，近似为0.67。如图7所示，虚线和实线箭头分别指示动脉和静脉血管。本文通过训练广义回归神经网络(GRNN)[10]分类器来分类动静脉血管，与其他分类器相比，其优势是不需要设定模型的表达式，直接通过核函数中的光滑因子调整网络性能，此外它学习快速，且当样本数量非常大的时候能够快速逼近。

实验选用DRIVE[11]视网膜图像库测试血管分类性能，此数据库有40幅彩色视网膜图像，该图库被分为测试集和训练集，每个集合中包含20幅图像。首先提取样本特征训练分类器，选训练集中的图像为训练样本，测试集中的图像作为测试样本。提取每一幅图像血管骨架上的点对应的11个特征(如表1所示)，将动脉骨架线用0标记(正样本)，将静脉骨架线用1标记(负样本)，在绿色通道图像上分别用黑色、白色显示。

利用训练好的GRNN分类器对测试集图像进行测试，DRIVE库中第2幅图像的分类结果如图8所示。动静脉血管分类为AVR计算打下基础。

图7 动静脉血管区分图

表1血管特征

特征标号特征 1骨架线上每个像素点归一化的管径值 2骨架线上每个像素点绿色通道归一化的灰度级 3骨架线上每个像素点红色通道归一化的灰度级 4-6骨架线上每个像素点绿色通道归一化的对比度、平均色调、饱和度 7骨架线上每个像素点绿色通道归一化的亮度 8-11标准偏差为2, 4, 6, 8的高斯函数对绿色通道图像进行滤波后骨架线上对应像素点的灰度级

6 AVR计算

医学上，视网膜血管管径测量的感兴趣区域是距离视盘边缘2~3个视盘半径内的圆环区域，之所以选择该区域是由于该区域的血管重叠少，可靠性高[1]，因此准确地定位视盘将有助于确定测量感兴趣区域并计算AVR。本文采用基于相位一致性的视网膜图像视盘定位方法[12]定位视盘，在感兴趣测量区域内的每一条血管划出5个等距离的测量点，如图9(a)所示。测量每段血管分支上5个等距测量点的管径并计算每段血管分支的管径平均值，静脉和动脉血管管径平均值分别记为和，如图9(b)所示。白色骨架上的测量点为动脉测量点，黑色骨架上的测量点为静脉测量点；接着以视盘中心为中心，以3倍为半径，沿逆时针方向寻找与每段静脉骨架线最近的动脉血管骨架，即找到动、静脉血管对，计算动静脉血管对的AVR。

7 实验结果分析

由于医生需要存在病变的眼底血管管径数据来辅助分析病变情况，为了验证本文方法的性能，采用REVIEW数据库和DRIVE数据库两种图像库对本文方法进行测试。而 REVIEW数据库中包含了眼底血管病变图像集，DRIVE数据库中包含了正常眼底血管图像，所以为了验证本文方法的性能和临床实用价值，本文选取含有眼底血管病变图像集的REVIEW库中公开的管径测量数据作为金标准来验证本文管径测量方法的稳定性和准确性，进而在含有正常眼底血管图像的DRIVE库中测试本文的血管分类方法及AVR测量的有效性。

图8 动静脉血管分类结果    图9 血管管径测量感兴趣区域

图10 眼底图像分割结果

表2不同方法对DRIVE数据库中血管分割的ACC, TPR和FPR(%)

方法ACCTPRFPR 文献[13]93.8271.202.76 文献[14]94.7773.542.11 本文方法95.8472.861.94

7.1 血管分割性能分析

由于本文是在二值血管的基础上测量管径，因此本文方法血管分割效果至关重要。为说明本文血管分割结果的有效性，对DRIVE视网膜库进行实验，部分实验结果如图10所示。以专家的手动标识结果为金标准，从准确度(ACCuracy, ACC)、正确分割率(True Positive Rate, TPR)以及错误分割率(False Positive Rate, FPR) 3方面评价本文方法的分割性能，本文血管分割方法为无监督的方法，与近几年的无监督的分割方法对比，结果如表2所示。可以看出本文方法ACC和FPR均优于其他方法，TPR仅略比文献[14]低0.68%，原因在于末梢极细小血管未能分割出。管径测量的区域不涉及末梢细小血管，整体上本文方法对血管的分割效果较好，性能优于其他方法，能够为测量血管管径提供基础。

7.2血管管径测量性能

首先为了说明管径测量的性能，本文从稳定性和准确性来分析。由于REVIEW数据库中有公开的眼底图像血管片段以及对应的管径测量数据[15]，为了对比说明，对此数据库提供的血管片段进行管径测量。REVIEW数据库中包含3个图像集，分别是中央光反射图集(Central Light Reflex Image Set, CLRIS)，有2幅早期动脉粥样硬化眼底图像，选取了21个呈现强烈的中央光反射的血管片段；分叉点图像集(Kick Point Image Set, KPIS)，有2幅图像，每幅图像中含有一段包括分叉的血管段；眼底血管病图像集(Vascular Disease Image Set, VDIS)，有8幅存在眼底血管病变的眼底图像，选取了79个存在病理噪声的血管片段。本文以KPIS图像集为例，部分结果如图11所示。此数据库包含3个观察者对不同形态的血管管径测量数据，以3个观察者数据的平均值作为金标准。

为验证本文方法管径测量结果的稳定性，本文对REVIEW数据库中的KPIS, CLRIS和VDIS图集进行测试。从管径平均值，管径测量平均值与金标准绝对误差和管径标准偏差两个方面将本文方法与其他方法进行对比，结果如表3所示。从表3中可以看出，本文方法得到的绝对值分别为0.04, 0.05, 0.01,分别为0.15, 3.74, 2.41，可以看出本文方法管径测量值可以得到较低的和，从而可以验证本文管径测量的稳定性。

图11 KPIS图集边界测量标记结果

7.3动静脉分类以及AVR性能

为验证动静脉血管分类的有效性，与眼科专家手动标识的DRIVE库测试集中的视网膜图像动静脉血管进行对比，结果如图12所示，图12(a)为眼科专家手动标识效果，白色标识为静脉，黑色标识为动脉。图12(b)为本文方法分类结果，可以看出整体上本文方法对动静脉血管能够准确分类，但对末梢血管段上的部分血管像素存在错分现象(图中白色方框标出)，原因在于这些部位的血管特征不明显导致错分。

表3不同方法对REVIEW库中管径测量数据对比(管径平均值以及不同方法与金标准的差，单位：像素；标准偏差)

方法KPIS数据对比CLRIS数据对比VDIS数据对比 金标准7.40 -0.3613.80 -4.128.85-2.57 文献[5]5.87-1.53- 7.00 -6.80-6.59-2.26- 文献[6]7.54 0.140.2413.8003.898.87 0.022.22 本文方法7.36-0.040.1513.75 -0.053.748.84-0.012.41

图12 专家手动分类和本文方法分类结果

将本文方法对DRIVE库中测试集中的20幅图像分类结果与专家手动标识结果对比，并绘制出ROC(Receiver Operating Characteristic)曲线，横轴为特异性即真正类率，代表静脉分类的准确率，纵轴为灵敏度即真负类率，代表动脉分类的准确率，如图13所示。其中实线和虚线分别为GRNN和K近邻分类对应的ROC曲线，GRNN分类ROC曲线下的面积为0.87, K近邻分类ROC曲线下的面积为0.84，从ROC曲线来看，相比K近邻，本文分类精度提高了3%。

为了验证ROI内AVR测量结果的有效性，统计了DRIVE训练集中图像的AVR测量结果。图14给出了本文AVR测量结果与眼科专家测量结果的对比图，横轴为本文方法动静脉血管平均比值(AVR)，纵轴为本文方法测量值与眼科专家的误差，上下两条线段为AVR误差的界限。

根据图14中的测量结果，通过计算95.35%的AVR值绝对误差在0.15以内，验证了AVR测量的有效性。

总之，本文方法对于DRIVE库进行血管分割的准确度和错误分割率均优于现有其他方法，正确分割率与其他方法相比也处于较高水平；对于REVIEW数据库进行管径测量的绝对值和均低于现有其他管径测量方法；对于DRIVE库进行动静脉血管分类和AVR测量，与眼科医生的手动标记做对比也具有较高的有效性。

8 结论

本文提出了一种基于Hesse矩阵和多尺度分析的动静脉血管管径自动测量方法，并设计了一套视网膜血管动静脉管径及AVR测量方案，包括血管分割、管径测量、动静脉血管分类和AVR测量。首先对DRIVE库进行测试，验证了本文血管分割的有效性；然后对REVIEW库中图像测试以及与其公开的管径测量数据对比，实验结果表明本文方法能够保证管径测量的稳定性和准确性，且能很好地作用于血管分叉和交叉处；最后提取血管管径等特征来训练GRNN分类器，对动静脉血管进行分类，进而在此基础上准确测量AVR。实验结果表明了本文方法的有效性。

图13 动、静脉血管分类ROC曲线          图14 AVR测量结果对比

[1] 李春艳, 徐力. 糖尿病视网膜血管管径的研究进展[J]. 医学综述, 2015, 25(17): 3173-3176. doi: 10.3969/j.issn.1006-2084. 2015.17.038.

LI Chunyan and XU Li. Research progress of diabetic retinal vascular caliber[J],, 2015, 25(17): 3173-3176. doi: 10.3969/j.issn.1006-2084.2015.17.038.

[2] NIEMEIJER M, XU X, DUMITRESCU A V,. Automated measurement of the arteriolar-to-venular width ratio in digital color fundus photographs[J]., 2011, 30(11): 1941-1950. doi: 10.1109/TMI.2011.2159619.

[3] GRISAN E and RUGGERI A. A divide et impera strategy for automatic classification of retinal vessels into arteries and veins[C]. The 25th Annual International Conference of the IEEE in Medicine and Biology Society, Cancun, Mexico, 2003, 1: 890-893. doi: 10.1109/IEMBS.2003.1279908.

[4] 姚畅, 陈后金. 一种基于先验知识的视网膜血管宽度测量方法[J]. 中国科技论文在线, 2009, 4(1): 64-68. doi: 10.3969/ j.issn.2095-2783.2009.01.011.

YAO Chang and CHEN Houjin. Measurement of retinal vessel widths based on prior knowledge[J]., 2009, 4(1): 64-68. doi: 10.3969/j.issn.2095-2783.2009. 01.011.

[5] LOWELL J, HUNTER A, STEEL D,. Quantitative measurement of retinal vessel widths from fundus images based on 2D modelling[J]., 2004, 23(10): 1196-1204. doi: 10.1109/TMI.2004. 830524.

[6] LUPA S CU C A, TEGOLO D, and TRUCCO E. Accurate estimation of retinal vessel width using bagged decision trees and an extended multiresolution Hermite model[J]., 2013, 17(8): 1164-1180. doi: 10.1016/j.media. 2013.07.006.

[7] FRANGI A F, NIESSEN W J, VINCKEN K L,. Multiscale Vessel Enhancement Filtering[M]. Berlin, Germany, Medical Image Computing and Computer-Assisted Interventation, Springer, 1998: 130-137.

[8] 张新红, 张帆, 崔延斌. 基于多尺度自适应滤波的DSA血管增强[J]. 计算机工程与应用, 2015, 51(14): 179-185. doi: 10. 3778/j.issn.1002-8331.1406-0457.

ZHANG Xinhong, ZHANG Fan, and CUI Yanbin. Blood vessel enhancement algorithm for DSA images based on adaptive multi-scale filtering[J]., 2015, 51(14): 179-185. doi: 10.3778/j.issn. 1002- 8331.1406-0457.

[9] CHEN Y, BARCELOS C A Z, and MAIR B A. Smoothing and edge detection by time-varying coupled nonlinear diffusion equations[J]., 2001, 82(2): 85-100. doi: 10.1006/cviu. 2001.0903.

[10] 严颂华, 吴世才, 吴雄斌. 基于神经网络的高频地波雷达目标到达角估计[J]. 电子与信息学报, 2008, 30(2): 339-342.

YAN Songhua, WU Shicai, and WU Xiongbin. DOA estimation based on neural network for HFGWR[J].&, 2008, 30(2): 339-342.

[11] STAAL J, ABRAMOFF M D, NIEMEIJER M,. DRIVE [OL]. http://www.isi.uu.nl/Research/Databases/DRIVE/, 2004.

[12] LI X, XIAO Z, and ZHANG F,. Fundus optic disc localization method based on phase congruency[J]., 2013, 7(4): 726-733.

[13] ZHANG B, ZHANG L, ZHANG L,. Retinal vessel extraction by matched filter with first-order derivative of Gaussian[J]., 2010, 40(4): 438-445. doi: 10.1016/j.compbiomed.2010.02.008.

[14] ZHAO Y Q, WANG X H, WANG X F,. Retinal vessels segmentation based on level set and region growing[J]., 2014, 47(7): 2437-2446. doi: 10.1016/j. patcog.2014.01.006.

[15] AL-DIRI B, HUNTER A, and STEEL D. Review[OL]. http://reviewdb.lincoln.ac.uk/, 2008.

Measurement of Retinal Diameters of Artery and Vein Based on Hesse Matrix and Multi-scale Analysis

XIAO Zhitao CUI Ning WU Jun GENG Lei ZHANG Fang WEN Jia TONG Jun LIU Xiaoting YANG Song

(,,300387,)

Many systemic diseases can cause changes of the diameters of retinal vessels and Arteriolar-to-Venular diameter Ratios (AVR), so it is of great importance to make quantitative analysis of the diameter of retinal vessels accurately in the diagnosis of the disease. An automatic method measuring the diameters of the artery and the vein and the AVR is proposed. Firstly, based on the segmentation of vascular network, the diameters of retinal vessels are measured according to the advantage of Hesse matrix for detecting line-like structure, and accurate localization of vascular direction with multi-scale analysis. Secondly, a General Regression Neural Network (GRNN) classifier is used to classify the artery and the vein points. Finally, the AVR in the Region Of Interest (ROI) is calculated. The validity of the proposed method is demonstrated by testing on the DRIVE and the REVIEW database.

Retinal vessels; Hesse matrix; Diameter measurement; General Regression Neural Network (GRNN); Arteriolar-to-Venular diameter Ratios (AVR)

TP391.4

A

1009-5896(2016)11-2871-08

10.11999/JEIT160165

2016-02-24；改回日期：2016-06-22；

2016-09-08

吴骏 zhenkongwujun@163.com

国家自然科学基金(61401439)，天津市科技支撑计划重点项目(14ZCZDGX00033)，天津市应用基础与前沿技术研究计划项目(15JCYBJC16600)，高等学校博士学科点专项科研基金(20131201110001)

The National Natural Science Foundation of China (61401439), Tianjin Science and Technology Supporting Projection (14ZCZDGX00033), Tianjin Research Program of Application Foundation and Advanced Technology (15JCYBJC 16600), Specialized Research Fund for the Doctoral Program of Higher Education of China (20131201110001)

肖志涛： 男，1971年，教授，研究方向为图像处理与模式识别、智能信号处理技术与系统.

崔 宁： 女，1990年，硕士，研究方向为图像处理与模式识别.

吴 骏： 男，1978年，副教授，研究方向为图像处理与模式识别、人工神经网络.