薛岚燕, 余 轮, 林嘉雯, 曹新容, 郑绍华

(1. 福州大学物理与信息工程学院， 福建 福州 350116； 2. 福建农林大学计算机与信息学院， 福建 福州 350002)

0 引言

视网膜血管可用非侵入方法直接观察， 为观察全身性血管异常的具体表现提供了客观条件， 使医生通过检查眼底发现或监测这类疾病成为可能， 因此， 眼底检查一直是全身血管性疾病简便、 经济的筛查与监测工具之一[1]. 视网膜循环系统与脑部循环系统有着相同的解剖生理结构， 而脑血管疾病是一类主要以小血管病变为基础的疾病[2]. 以高血压和脑卒中等疾病为例， 眼底观察表现为视网膜微血管改变， 高血压主要与局限性动脉缩窄、 动静脉交叉压痕有关[3]， 脑卒中主要与视网膜动脉缩窄、 视网膜静脉扩张等特征有关[4]. 因此， 观察视网膜血管的改变为脑血管疾病的早期诊断、 甚至预测其危险性提供了理论依据[5].

图1 本算法过程图Fig.1 Framework of the proposed method

视网膜血管直径信息对于临床诊断具有重要意义， 很多科研人员提出了计算机自动管径测量算法[6]. 测量方法主要分为以下两种. 1) 基于血管轮廓的管径测量. Delibasis等[7]采用圆形结构估计法， 通过计算血管两侧在圆上的点数得到血管与圆心的夹角， 根据几何关系得到血管直径. Abdolhossein等[8]在提取血管网络和中心线基础上， 以中心线为圆心， 以经验值R为半径做圆， 通过半径和圆与血管边缘交点间的夹角计算血管直径. Kumar等[9]利用LDA分类器将血管横截面上的血管和背景进行分类， 并将横截面上的血管宽度作为血管管径. 林土胜等[10]结合邻域定位技术来自动确定边界点， 通过建立管径增量符号比较法来自动寻优求取管径. 基于血管轮廓的管径测量方法需要确切的血管网络边界识别， 因此对血管边界的准确性要求很高. 如果血管边界识别不准确会导致宽度测量误差. 2) 基于血管横截面灰度变化曲线的模型拟合测量. 国内外学者多数是采用高斯曲线来拟合血管横截面轮廓， 通过计算轮廓曲线半高度等宽来得到管径. Zhou等[11]提出了一种基于高斯曲线模型的视网膜血管直径测量方法. Carmen等[12]采用决策树和Hermite模型对视网膜血管宽度进行测量. 吴辉群[6]采用基于最小二乘法的混合高斯模型对血管横截面像素灰度值分布曲线进行拟合， 以半高度全宽计算方法获得血管直径. 基于血管横截面灰度变化曲线的模型拟合测量需要借助血管中心线做垂线获取血管横截面的曲线分布， 血管宽度测量与准确的血管网络、 所选取的拟合曲线以及所选定的血管断面灰度曲线位置均有密切关系. 本研究针对这两种方法存在的问题， 提出了一种交互式管径测量方法， 目的是通过交互式方法捕捉边界点， 克服边界选取的随机性， 并且对一个小血管段测量多个管径求平均值作为该血管段的管径值， 从而提高测量精度.

采用自适应对比度增强(ACE)对眼底图像进行增强， 引入盒式滤波用于ACE的计算， 不仅增强效果好， 同时降低了计算复杂度， 速度也得到了极大的提高. 通过2D Gabor滤波提取眼底血管轮廓， 进行阈值分割并细化提取骨架， 确定骨架上的分叉点和终点， 再对分叉点进行连通域判断确定最终血管网络. 最后采用交互式方法得到血管直径， 系统流程如图1所示.

1 眼底血管网络的提取

1.1 基于盒式滤波自适应对比度增强的眼底图像增强

眼底图像具有背景不均及血管对比度弱的特点， 为增强眼底血管改变的细微变化， 同时考虑到大批量眼底普查时机器自动识别的需要， 提高诊断精度， 对快速自动的眼底血管图像增强方法提出了迫切要求[13]. 采用自适应对比度增强(ACE)， 该算法充分考虑到局部信息并且可以限制噪声的放大， 具有很好的增强效果， 但存在计算量大的问题. 考虑到采用ACE对眼底图像进行增强， 图像的局部矩形内像素的和、 平方和、 均值、 方差等特征会频繁地在算法中使用， 因此需要对这些特征进行优化. 本算法采用盒式滤波(boxfilter)优化方法， 盒式滤波能做到与窗口大小无关， 使复杂度为O(MN)的求和， 求方差运算降低到O(1)或近似于O(1)的复杂度， 大大降低了计算复杂度并提高计算速度， 是常规均值滤波的一种有效的滑动平均求和加速算法.

为了通过增强信号以便于提取眼底血管， 需要对眼底图像中的高频成分进行放大. 在得到低频成分后， 将原图减去低频成分得到高频成分并对高频成分增加一个对比度增益从而放大高频成分， 并与低频成分相加得到增强的图像. 自适应对比度增强算法[14]描述如下， 假设x(i,j)为眼底图像中的某一点， 则图像的低频部分， 即局部的平均值可表示为：

(1)

则通过对高频成分放大得到增强的图像：

f(i,j)=mx(i,j)+G(i,j)[x(i,j)-mx(i,j)]

(2)

其中：G(i,j)为对比度增益放大系数， 可表示为：

(3)

其中：σx(i，j)为采用均值滤波器滤波前后的误差;D为图像的全局平均值.

其中：σ2x(i,j)为采用均值滤波器滤波前后的均方误差.

通过式(2)可以看出， 在整个算法中频繁出现图像的局部矩形内像素的和、 平方和、 均值、 方差等特征的计算， 计算量非常大， 大大影响了计算速度. 考虑到计算复杂度和速度， 本研究采用盒式滤波优化算法， 可以简化特征的计算. 它的基本原理为： 首先建立一个数组A， 宽高与原图像相等， 然后对这个数组进行初始化赋值， 每个元素的值A[i]赋为该像素邻域内的像素和(或像素平方和)， 在求某个矩形内像素和的时候， 直接访问数组中对应的位置即可. 具体过程如文[15]所描述.

采用自适应对比度增强并结合了盒式滤波优化算法进行运算， 具体过程描述如下：

1) 首先对原图img进行半径为n的盒式滤波运算. 由于n值过大过小都会造成图像负优化， 根据经验值选择n=50, 得到每个像素都是原图对应像素的邻域内的均值， 记为matrix_mean， 即matrix_mean=boxfilter(source=img,r=n)， 此步骤计算完所有像素点对应邻域内的均值.

古代中国的节日生活优雅、端庄，每个节日背后都有层叠累积、意蕴深长的故事传说。 在节日与节日的岁月轮回中，生命的异彩照亮了平淡的日子，焕发出生活的艺术和艺术的生活。

2) 计算原图img与matrix_mean的均方误差， 记为matrix_sigmas_total， 即matrix_sigmas_total=(img-matrix_mean)^2.

3) 对matrix_sigmas_total进行半径为n的boxfilter运算， 得到每个像素都是原图对应像素邻域内的方差， 开方后即为均方差(标准差)， 记为matrix_sigmas， 即matrix_sigmas=sqrt(boxfilter(source_sigmas_total,r=n)) .

4) 最后计算全局平均值D， 根据公式(2)计算出增强的图像灰度值. 通过将盒式滤波应用到自适应对比度增强算法中的计算， 由于盒式滤波的初始化过程非常快速， 每个矩形的计算基本上只需要一加一减两次运算， 在具体求某个矩形特征时， 会降低到1次加减运算. 图2对应的是原图和对比度增强图， 相对于原图， 经过增强后的图像中血管和背景对比度提高了， 有利于血管分割.

图2 自适应对比度增强对比图Fig.2 Adaptive contrast enhancing images

1.2 2D Gabor滤波

2D Gabor小波在进行视网膜血管位置特征检测时具有方向选择性， 而且对于不同宽度的视网膜血管， 2D Gabor小波可以选择不同尺度参数对其进行处理[16]. 因此, 采用2D Gabor滤波器作为核函数来提取眼底血管轮廓特征. 二维Gabor滤波函数可以表示为：

(4)

x0=xcosθ+ysinθ

(5)

y0=-xsinθ+ycosθ

(6)

图3 Gabor变换图Fig.3 Gabor image

1.3 基于骨架的血管网络提取

血管的宽度可以细至1 px， 宽至20 px(取决于血管的宽度和图像分辨率)， 即血管粗细的比例可达20倍[17]. 在对眼底图像增强和滤波处理后需要进行阈值分割得到二值图像， 可以消除图像灰度插值影响， 使血管和背景之间能形成明显边界以解决管径边界的确切性并消除量化的随机误差[18]. 由于一维阈值分割算法稳定性较差并且对噪声敏感， 而二维阈值分割算法计算复杂度高， 为了解决上述问题， 采用基于区域的改进一维直方图阈值分割算法[19]. 假设f(m,n)表示像素值，g(m,n) 表示其邻域，z为灰度级,z=0, 1, 2, …,L-1， 根据f(m,n)和g(m,n)的关系采用最大值统计方法计算直方图， 可以表示为：

(7)

(8)

以新直方图所得阈值对Gabor响应图进行分割. 由于Gabor滤波器特性， 提取结果中含有大量的非血管结构， 需要将二值分割后的血管细化为单像素宽高的血管候选骨架结构， 通过确定骨架上血管的分叉点和终点来跟踪抑制非血管结构从而确定真正的血管结构[20]. 血管在图像中呈现连通的网络结构， 长度小于10 px的分支属于其他非血管区域残留的可能性很大， 可以考虑去除. 具体步骤如下：

1) 通过直方图阈值法对Gabor响应图进行二值分割.

2) 对二值图像进行细化处理得到单像素宽的血管候选骨架， 其中细化处理采用Hilditch细化算法[21]来实现.

3) 对细化后的图像标记分叉点和终点， 遍历图像上所有像素值为1的像素点. 判断其八邻域的连通数， 如果连通数为1， 则判断为终点； 如果连通数等于2， 则判断为普通点； 如果连通数大于2， 则判断为血管交叉点. 图4为血管网络提取过程示意图.

图4 血管网络提取过程示意图Fig.4 Extracting process of the vessel network

4) 在细化图像上随机选取一个交叉点， 采用种子填充法判断其连通区域并进行区域填充. 假如没有遇到断点或终点， 但是超出了设定的缓存数组极限， 则停止； 假如没有超出缓存， 但是遇到了终点， 那么表示这段血管填充结束； 如果在得到的连通域出现过多分叉(通常只有很细小的血管才会有过多分叉)， 则丢弃超出分叉数的部分.

5) 重复4)的步骤， 直到所有交叉点判断完成.

2 交互式最优管径测量法

视网膜血管直径测量方法是利用数字化视网膜图像， 以交互式方法来测量感兴趣区域内所有可测量的小动脉和小静脉， 其中感兴趣的测量区域[22]是以视盘中心为圆心， 1 DD ( DD表示视盘直径)和1.5 DD分别画圆， 两圆中间部分定义为感兴趣测量区域. Knudtson等[23]进一步将测量方法细化至 6 条最大的小动脉和 6 条最大的小静脉来测量视网膜血管， 每条血管总共测量 5 次， 选定血管段每次测量间距尽量相等， 一条血管得到5个子血管段. 将五个子血管段得到的管径求平均即为该条血管上的管径. 因此， 在获得血管网络后采用交互式方法对感兴趣区域的动静脉测量管径[24-25].

1) 在血管两侧采集4个点， 采集顺序采用N字序列采集， 即血管两侧的两点顺序一致即可(相对于U字序列， 两边采集顺序相反). 根据血管两侧的两点分别自动计算两点间的1/3和2/3位置处的点， 从而两边共同扩展成8点， 之所以只通过手动点击4个点而不是8个点是防止可能不满足中垂线必须穿过对面两点之间的要求， 如图5所示， 在血管两边按照某种顺序依次选取A′至H′ 8个点. 选取顺序决定了后面血管另一侧对应点的选择， 则A′对应E′.

2) 利用这8个点来寻找边界， 将血管两侧同一水平线上的两点进行连线， 这两点分别沿着连线方向向血管边界方向移动并计算经过的每个像素点的像素值. 当像素值非 0 时， 则认为到达边界像素点. 通过该步处理使8个血管外的点变为8个在血管边界上的点. 如图6所示， 将A′ 移至A， 将E′ 移至E， 以此类推.

图5 序列采集图Fig.5 Sequence image

图6 寻找边界点图Fig.6 Searching edge point

图7 测量管径图Fig.7 Measuring the diameter

3) 计算紧邻两点间的中点到对面对应两点线段与该中点的中垂线的交点距离记为一次计算. 如图7所示， 取线段AB， 还有线段AB中点X， 过X作AB的垂线， 与线段EF交叉于点X′(注意不是中点), 则XX′的长度为一个管径长度. 图7中有A、B、C、D、E、F、G、H等8个点， 以此类推共有6个XX′这样的管径长度， 取其平均值作为最后的测量值.

在上述步骤中需要注意， 在计算交叉像素点与过某一侧血管两点的距离的平均值时需要加上修正值. 这是因为在选取的边缘点位于血管内部， 计算两点间距离得到的并不是血管的管径. 通过实验发现, 在测量垂直或水平的血管少了1 px， 45°方向的少了1.40 px， 30°方向的少了1.15 px. 因此, 在设置测量垂直或水平的血管时修正值设为1 px， 其他为1.20 px.

3 实验结果与分析

本算法是基于Visual Studio2008和OpenCV2.4.3， 在主频可达3.8 GHz， CPU为intele31231V3， 内存16 GB的计算机上进行实验， 编程语言使用C++. 分别采用DRIVE和STARE眼底图库以及医院临床图像共150幅眼底图像进行实验和性能评估.

3.1 盒式滤波自适应对比度增强结果

一幅大小为MN的眼底图像， 采用自适应对比度增强优化算法进行处理， 相对于传统算法， 在增强效果上保持一致， 但由于盒式滤波算法中每个矩形的计算基本上只需要一加一减两次运算， 在具体求某个矩形特征时， 会降低到1次加减运算， 因此大大降低了计算复杂度和提高了计算速度. 不同尺寸的眼底图采用两种算法的计算时间进行比较， 采用不同大小的窗口半径计算所需要的时间， 结果如表1所示. 从表1中数据可以得出， 对同一个尺寸的眼底图像， 随着窗口半径的增加， 传统算法的计算时间数倍增加， 但优化算法的计算时间变化不大， 在计算速度上要比传统算法提高很多.

表1 不同尺寸眼底图像两种算法计算时间比较结果

3.2 眼底血管网络提取的结果

为验证血管分割算法， 分别计算TPR、 FPR和ACC:

其中: TP和TN分别表示正确的血管像素和正确的背景像素， 即与眼科专家判断一致; FP和FN分别表示错误的血管像素和错误的背景像素， 即与眼科专家判断不一致; TPR为灵敏度， 表示正确分割的血管像素点数占血管像素点总数的比值; FPR为特异性， 表示错误的血管像素点数占背景像素点总数的比值. ACC为准确率， 表示正确分割的像素点数占像素点总数的比值.

不同图像的血管网络提取结果如图8所示， 第一列代表眼底图像原图， 第二列代表血管填充图, 第三列代表血管网络图. 在对眼底彩色图像增强和滤波处理后进行阈值分割得到二值图像. 由于提取结果中含有大量的非血管结构， 需要对二值分割后的血管进行细化为单像素宽高的血管候选骨架结构， 通过确定骨架上血管的分叉点和终点进行跟踪抑制非血管结构从而确定真正的血管结构.

图8 血管网络提取结果Fig.8 Extracting process of the vessel network

方法TPRFPRACC文[16]0.75260.02570.9457文[26]0.68980.03040.9416文[27]0.65200.02550.9461文[28]--0.9461本算法0.76120.02510.9475

表2给出了本算法和其他算法对DRIVE眼底图库进行实验所得到的灵敏度、 特异性和准确率的比较结果， 由表中数据可知， 本算法的准确率要高于其他算法， 可达94.75%。

3.3 血管管径测量结果

本研究共测试STARE数据集中 81 幅眼底图像， 对每幅图像感兴趣测量区域ROI中最粗的6条动脉和静脉， 采用交互式管径测量法测量管径， 如图9所示. 该图像为STARE库中im 0163.ppm眼底图像， 选取其感兴趣测量区域(图像中的坐标系)中某一段血管， 如图9(b)上黄框所示， 将该区域放大， 图9(c)至图9(g)为交互式管径测量过程示意图. 对感兴趣测量区域中的6条最粗的动静脉(如图9(a)所示， 蓝色粗框为静脉， 红色细框为动脉)分别通过交互式管径测量法、 人工确定血管直径法以及文[6]中所提到的基于血管横截面灰度曲线混合高斯曲线拟合的宽度测量算法进行测量[20]， 三种方法测量的血管直径分别记为D0,D1,D2, 其中人工确定血管直径法是眼科医生直接对眼底图像中的血管边界进行标识再测量血管的直径， 测量结果如表3所示.

从表3中数据可见， 本算法与另外两种方法进行比较， 三种方法所测量的血管直径基本一致， 误差很小； 对整个数据集感兴趣测量区域中的血管段进行管径测量， 与人工确定血管直径法的平均误差为0.032 px， 与文 [6]中的算法的平均误差为0.045 px， 两组算法结果有较高的一致性.

图9 交互式最优管径测量过程示意图Fig.9 Measurement process of the interaction optimal diameter

4 结语

本研究在利用结合盒式滤波自适应对比度增强优化算法提高图像质量的基础上， 提取血管网络并利用交互式测量方法对血管直径进行测量. 实验表明， 血管分割方法在灵敏度、 特异性和准确率方面较其他分割方法有所提高， 在感兴趣测量区域内实现血管直径测量， 测量结果与人工测量结果接近， 准确性高. 因此， 方法对于视网膜血管改变的早期诊断具有一定的实用价值.

参考文献：

[1] 王爽, 徐亮, 李建军. 视网膜微血管异常与心脑血管疾病关系的流行病学研究[J]. 国外医学(眼科学分册), 2005, 29(3): 145-148.

[2] 汪晓磊, 满凤媛, 王艳玲， 等. 基于视网膜血管评估心脑血管疾病的研究进展[J]. 临床和实验医学杂志, 2015, 14(22): 1921-1923.

[3] 刘雪, 王爽, 徐亮. 视网膜血管改变与高血压的相关性[J]. 国际眼科纵览, 2015, 39(1): 1-7.

[4] 张莉, 徐亮, 杨桦， 等. 眼底指标改变与脑卒中患病的相关性[J]. 眼科, 2015, 24(1): 13-18.

[5] PATTON N, ASLAM T, MACGILLIVRAY T,etal. Retinal vascular image analysis as a potential screening tool for cerebrovascular disease: a rationale based on homology between cerebral and retinal microvasculatures[J]. J Anat, 2005, 206(4): 319-348.

[6] 吴辉群. 慢性病信息管理系统中视网膜图像的互操作性及其血管网络定量分析研究[D]. 上海： 复旦大学, 2014.

[7] DELIBASIS K K, KECHRINIOTIS A I, TSONOS C,etal. Automatic model based tracing algorithm for vessel segmentation and diameter estimation[J]. Computer Method Programs in Biomedicine, 2010, 100(2): 108-122.

[8] ABDOLHOSSEIN F, AHMAD R N. Automatic wavelet-based retinal blood vessels segmentation and vessel diameter estimation[J]. Biomedical Signal Processing and Control， 2013(8): 71-80.

[9] KUMAR D K, ALIAHMAD B, HAO H. Retinal vessel diameter measurement using unsupervised linear discriminate analysis[J]. ISRN Ophthalmol, 2012, 2012: 1-7.

[10] 林土胜, 徐锦堂. 测量视网膜血管管径的新方法[J]. 眼科研究, 1998, 16(3): 231-232.

[11] ZHOU L, RZES Z O, TARSKI M S,etal. The detection and quantification of retinopathy using digital angiograms[J]. IEEE Transactions on Medical Imaging, 1994, 13(4): 619-626.

[12] LUPASCU C A, TEGOLO D, TRUCCO E. Accurate estimation of retinal vessel width using bagged decision trees and an extended multiresolution Hermite mode[J]. Med Image Anal， 2013, 17(8): 1164-1180.

[13] 许雷, 郑筱祥. 一种基于小波变换及数学形态学方法的眼底图像增强即定量分析方法[J]. 生物物理学报, 1998, 14(1): 70-76.

[14] SANTHI K, BANU R S D W. Adaptive contrast enhancement using modified histogram equalization[J]. International Journal for Light and Electron Optics, 2015,126(19): 1809-1814.

[15] MCDONNELL M J. Box-filtering techniques[J]. Computer Graphics & Image Processing, 1981,17(1): 65-70.

[16] 王晓红, 赵于前, 廖苗, 等. 基于多尺度2D Gabor小波的视网膜血管自动分割[J]. 自动化学报, 2015, 41(5): 970-980.

[17] FRAZ M M, REMAGNINO P, HOPPE A,etal. Blood vessel segmentation methodologies in retinal images: a survey[J]. Computer Methods and Programs in Biomedicine, 2012, 108(1): 407-433.

[18] 林土胜, 徐锦堂, 梁庆源. 视网膜血管管径量化及对眼底血管病变的应用价值[J]. 暨南大学学报(医学版), 1999, 20(2): 33-36.

[19] 汪启伟. 图像直方图特征及其应用研究[D]. 合肥： 中国科学技术大学, 2014.

[20] 李丽华, 王凯. 基于数学形态学的视网膜血管提取算法[J]. 北京生物医学工程, 2014, 33(5): 497-501.

[21] YU J, LI Y. Improving hilditch thinning algorithms for text image[C]//International Conference on E-Learning, E-Business, Enterprise Information Systems, and E-Government. HongKong: IEEE Computer Society, 2009: 76-79.

[22] NIEMEIJER M, XU X Y, DUMITRESCU A V,etal. Automated measurement of the arteriolar-to-venular width ratio in digital color fundus photographs[J]. IEEE Transactions on Medical Maging, 2011, 30(11): 1941-1950.

[23] KNUDTSON M D, LEE K E, HUBBARD L D,etal. Revised formulas for summarizing retinal vessel diameters[J]. Current Eye Research, 2003, 27(2): 143-149.

[24] FATHI A, NAGHSH-NILCHI A R. Automatic wavelet-based retinal blood vessels segmentation and vessel diameter estimation[J]. Biomedical Signal Processing & Control, 2013, 8(1): 71-80.

[25] 林土胜, 张莺, 谢洁珍. 视网膜血管网络的边界识别和寻优算法[J]. 电路与系统学报， 1997, 2(4): 62-64.

[26] NIEMEIJER M, STAAL J, GINNEKEN B V,etal. Comparative study of retinal vessel segmentation methods on a new publicly available database[C]//Proceedings of SPIE-The international society for optical engineering. Canberra: IEEE, 2004, 5370: 648-656.

[27] KABA D, SALAZAR-GONZALEZ A G, LI Y,etal. Segmentation of retinal blood vessels using gaussian mixture models and expectation maximisation[M]. Heidelberg: Springer, 2013: 105-112.

[28] WANG Y F, JI G R, LIN P,etal. Retinal vessel segmentation using multiwavelet kernels and multiscale hierarchical decomposition[J]. Pattern Recognition, 2013, 46(3): 703-715.