1 引言
对视网膜血管的参数进行准确估计一直是眼底图像分析中重点关注的问题. 因为大量临床经验表明，视网膜的动脉和静脉宽度变化与糖尿病、 心脑血管疾病等都密切相关[1,2]，如动脉变窄与高血压有关； 静脉变窄则患有中风和冠状动脉粥样硬化性心脏病的可能性较大. 动静脉宽度之比(arteriolar-to-venular diameter ratio，AVR)可用来衡量视网膜血管的宽度变化，但其相对变化十分微小，在临床检查时即使是对于经验丰富的眼科医生也是十分困难的，即从眼底图像中人工或半自动地估计AVR值是复杂而费力的. 因此，利用图像处理和计算机辅助方法从眼底图像中自动地测量AVR具有十分重要的临床意义，有助于医生通过眼底图像对病人的患病危险进行评估，实现对中风、 认知能力下降、 脑萎缩、 心肌梗塞、 糖尿病视网膜病变、 早产儿视网膜病变等疾病的筛查、 预防和早期治疗.
国内外许多学者都对眼底图像中AVR的计算展开了深入研究. Hubbard等人[3, 4]提出了动静脉直径的半自动计算法，并将其应用于社区动脉粥样硬化风险研究(atherosclerosis risk in communities，ARIC)，该方法首先将血管分为动、 静脉，然后在测量区域对动静脉进行分段和相应的直径计算； Patker等[5]提出了半自动计算方法，其输入为人工辨识出的血管对，通过自动检测血管壁的边缘计算出血管直径和AVR； Tramontan等[6, 7]提出了带有手动校正选项的AVR自动估计法，该方法首先通过血管跟踪确定了血管网络，之后在血管网络上定位了视盘，通过中心反射的特征对动静脉进行分类，使用直径大于45 μm的血管进行了AVR计算，其中血管的分割和分类可以通过手工进行校正； Nam等[8]的AVR计算方法是基于圆的光强分布提出的，该方法首先使用RGB(red green blue)图像中的红色通道进行视盘检测，然后以视盘中心为圆心，以视盘直径为半径作一个圆，沿着该圆可以得到绿色通道的强度分布，该分布的凹陷部分可视为血管，并根据形状对其进行动静脉的分类，血管宽度大于7个像素的用于计算AVR； 以上方法均是半自动方法，需要眼科医生的参与，其经验和技术对估计结果的影响较大. Niemeijer、 Muramatsu等[9, 10]提出了全自动AVR计算方法，该方法依次进行视盘的定位、 感兴趣区域(region of interest，ROI)的确定、 动静脉血管的分类、 AVR的计算，本文也采用了类似的过程，但每一步所使用的方法均不相同，尤其是在动静脉分类和血管宽度计算部分，本文提出的方法使得估计结果更加理想. 通过本文方法与其它算法在血管分割、 动静脉分类、 AVR测量等方面的性能对比，验证了本文算法的有效性、 准确性和一致性.
2 本文方法准确的血管宽度估计需要血管分割以分离血管与背景，需要检测视盘以确定进行计算和检测的感兴趣区域，需要准确的动静脉分类以测量血管宽度. 同时由于在眼底图像中视网膜血管比背景黑，因此首先要对采集到的图像进行预处理，对图像进行增强，以提高图像中的信噪比，使图像清晰化，减轻后续处理的压力. 因此，本文方法的具体过程如图 1所示.
 |
| 图 1 本文方法的计算过程Fig. 1 Measurement process of the proposed approach |
2.1 预处理
眼底图像在采集过程中往往受到非理想成像条件限制而使采集到的图像对比度不够强，严重影响参数自动检测，为此必须首先进行预处理以满足眼底血管病变检测系统的要求.
由于血管具有不同的宽度和方向，同时考虑到血管具有线性结构且横截面可以用高斯函数建模，所以这里采用连续小波变换(continuous wavelet transform，CWT)进行多尺度分析[11]. 在小波变换中，粗血管可以由低分辨率的滤波器得到，细血管可由高分辨率的滤波器得到，同时CWT的方向选择性可以用于提取不同方向上的血管. 为克服诸如亮斑等非血管结构的边缘对高斯核函数的强响应，本文选取了复Morlet母小波：
要将血管进行动静脉的区分首先要从眼底图像中将血管分割出来. 本文在眼底图像的绿色通道中采用Top-hat变换和双环滤波技术[12, 13]实现血管分割. 在对图像进行Top-hat变换时选取平面圆盘形结构元素，根据对训练集的测试，结构元素的尺度选为19个像素比较适合； 在进行双环滤波时，双环滤波器的内环半径选为7个像素，外环半径选为17个像素，使用这样尺度固定的双环滤波器对图像不同位置进行卷积，之后将Top-hat变换和双环滤波的结果相加，并通过设定置信度提取出血管(这里的置信度设置为15％). 通过形态结构去除了无需参与AVR计算的小于5个像素的细血管.
2.2.2 感兴趣区域的确定为确定感兴趣区域，必须要对视盘(optic disk，OD)进行定位，这里采用的是文[14]的方法. 视盘是眼底图像中较亮的一部分，可以近似为圆形. 因此该方法首先采用形态学处理方法进行4次开、闭迭代运算获得图像中最亮部分； 然后通过两步自动阈值选取方法找到OD及其附近的像素并去除OD附近较黑的像素，这两个阈值的选取都采用Otsu方法[15]； 最后，应用CHT变换(circular hough transform)获得OD的边界，从而得到了OD的中心点. 以OD的中心点为圆心、 1DD(disk diameter)和1.5DD分别画圆，两圆中间部分定义为本文的ROI(图 2中的B区).
 |
| 图 2 OD及ROI的确定Fig. 2 Definition of OD and ROI |
2.2.3 动静脉分类
动静脉血管分类的方法有自动和半自动两种，本文采用的是基于对主要血管进行特征识别的自动方法. 该方法由两个步骤完成： 第1步是选取适当的特征集，第2步是选择适当的分类方法，基于所提取的特征将分割出来的血管分为静脉和动脉.
在特征集选取之前，首先采用细化方法[16]分析提取出已分割的血管网络的中心线. 同时断开交叉或分叉的血管，具体做法是： 如果某个中心线像素的相邻像素个数多于两个则认为在这个中心线像素周围有交叉或分叉，因此在这些相邻像素方向寻找出这个中心像素点到非血管像素点的最短路径，将这些路径末端连接起来，去掉连接图形内部的像素后，交叉点及分叉点便被清除掉，继而达到断开交叉或分叉的血管的目的. 由于静脉是暗红色的，而动脉是鲜红色的，所以通过红色通道来区分动静脉是最佳选择. 但有些图像的红色通道已经非常饱和，提升的空间非常小，因此添加了绿色通道图像. 又由于一些图像绿色通道中动静脉的对比度非常低，我们还引入了LAB色彩空间和HSL空间. 为减少噪声的影响，以中心线上每个像素点为圆心做直径为10个像素的圆，在该区域计算灰度的均值、 方差和中值. 具体地，用下列等式表示特征向量，Vij表示基于中心线的特征向量.
在选择分类方法时，本文选择了无需监督学习的k均值聚类方法，其基本思想是通过迭代的方法，逐次更新各聚类中心的值，直至聚类中心不再发生变化，则说明数据对象全部分配到自己所在的聚类中，且聚类准则函数收敛(这里的聚类准则函数选为欧氏距离).
将上述血管特征向量作为算法的输入，则中心线的每一个像素都被划分为静脉或动脉中某一类，然后根据每一根血管中的动、 静脉像素的概率大小确定其类型，即计算式(9)、 式(10)：
 |
| 图 3 动静脉分类的结果Fig. 3 Results of artery/vein classification |
2.3 AVR计算 2.3.1 血管宽度测量
由于血管粗细是变化的，因此在计算宽度时，首先将血管进行分段，即从1DD到1.5DD以0.1DD为步长进行划分出6个圆环，从而ROI每个血管被分为5个血管子段(如图 4所示)，然后在每个血管子段的中心线上进行宽度测量. 计算宽度时，本文采用的是圆形结构估计法[17]，即以所分割出来的血管子段的中心点为圆心做一个半径为R(R=Wm+5，Wm是最大的血管宽度)的圆，通过计算血管两侧在该圆上的点数得到血管与圆心的夹角，继而由几何关系计算出血管宽度.
 |
| 图 4 血管子段的划分方法Fig. 4 Segments method of sub vessels |
2.3.2 计算AVR
AVR的计算公式是： AVR=CRAE/CRVE，其中CRAE(central retinal arterial equivalent，视网膜中央动脉直径)是选择ROI中6条最宽的动脉，应用Knudtson公式[18]计算得到的： CRAE=0.88×(wi2+wj2)，(i，j)是6条所选动脉中的最宽和最窄的动脉宽度； 类似地，CRVE(central retinal vein equivalent，视网膜中央静脉直径)是选择ROI中6条最宽的静脉进行计算： CRVE=0.95×(wi2+wj2)，(i，j)是6条所选静脉中的最宽和最窄的静脉宽度.
图 5给出了本文对CRAE和CRVE的迭代计算过程，向量A表示ROI中最宽的6条动脉的宽度，向量V表示ROI中最宽的6条静脉的宽度.
 |
| 图 5 CRAE和CRVE的计算流程图Fig. 5 Calculating flow chart of CRAE and CRVE |
按照上述对血管子段的划分方法，将会得到5个AVR值，这5个值的平均值为最终的AVR估计值.
3 实验结果 3.1 实验所用的图像(1) 何氏眼科数据库. 辽宁何氏眼科医院于2004年至2014年对15位(其中3位后来无法联系上)志愿者每年进行一次图像采集用于视网膜病变筛查系统的实验. 所采用的相机是YZ50A免散瞳眼底照相机，散瞳拍摄角度30°，小瞳孔拍摄角度45°，免散瞳拍摄角度27°，曝光量控制在30 ws～90 ws. 拍摄后得到1 400×1 225 pixel大小的JPEG格式眼部图像150幅，该筛查系统从中任选50幅图像，由权威医生对主要动静脉血管、 视盘、 视盘中心点分别进行了标注，可作为衡量标准. 其中30幅图像作为训练集，用于动静脉的分类和算法参数的确定，其余20幅作为测试集用于对算法和结果进行评估. 本实验中所有参数的设置均是基于此数据库.
(2) DRIVE数据库. 该数据库[19]是由荷兰的糖尿病视网膜病变筛查系统得到的，由40幅进行人工标注和分类的JPEG格式的视网膜图像组成，20幅作为训练图像，另外20幅作为测试图像. 图像所用相机为佳能的CRR5，图像大小是565×584 pixel. 第一位专家所标注的图像作为衡量标准.
3.2 血管分割结果本文的方法是基于像素的方法，每个像素被分类为血管或背景，因此对分类结果可产生4个参数： TP(true positive)表示像素在专家标注图像和本文分割算法中均被认定为血管； TN(true negative)表示像素在专家标注图像和本文分割算法中均被认定为非血管； FP(false positive)表示像素在专家标注图像中被认定非血管，而本文分割算法却将其认定为血管； FN(false negative)表示像素在专家标注图像中被认定血管，而本文分割算法却将其认定为非血管. 这4个参数可用来评价分割算法的好坏，如表 1所示. 表 2所示为本文算法对两个数据库中的20幅测试图像所得到的性能指标.
| 指标 | 描述 |
| Sensitivity | TN/(TN+FP) |
| Specificity | TP/(TP+FN) |
| Accuracy | (TP+TN)/(总的像素点数) |
| 指标 | 数据库 | Min | Max | Average |
| Specificity | 何氏 | 0.947 1 | 0.969 2 | 0.955 9 |
| Sensitivity | 何氏 | 0.648 5 | 0.847 6 | 0.756 8 |
| Accuracy | 何氏 | 0.920 1 | 0.948 5 | 0.938 6 |
| Specificity | DRIVE | 0.925 7 | 0.950 1 | 0.932 4 |
| Sensitivity | DRIVE | 0.638 7 | 0.839 8 | 0.732 5 |
| Accuracy | DRIVE | 0.895 4 | 0.921 6 | 0.918 4 |
Specificity表征了算法对非血管像素的检测能力，Sensitivity表征了算法对血管像素检测的能力，Accuracy表征了算法正确分类的能力. 图 6给出了对训练集进行血管分割时，Top-hat滤波器中结构元素的不同尺度选取对TP和FP的影响，可以看到由于本文所采集的眼底图像中最大的血管宽度大概为15个像素，因此Top-hat滤波器的结构元素尺度为11个像素时，所得的性能曲线处于较低位置，最好的结果是结构元素尺度取为19个像素.
 |
| 图 6 Top-hat滤波器中不同的结构元素尺度值对血管分割的结果影响Fig. 6 he result influence of the Top-hat filters with different diameters on vessel segmentation |
3.3 动、 静脉分类结果
在动、 静脉分类中，Sensitivity表示可以正确地将血管分为动脉或静脉的概率，动脉的Specificity表示对静脉进行正确区分的概率，静脉的Specificity表示对动脉进行正确区分的概率. 表 3对比了本文方法与现有两种方法的分类结果，可以看到本文静脉的Sensitivity和Specificity(分别用Sensv和Specv表示)与文[20]和文[21]的结果基本相当，且优于文[21]的方法，而本文方法得到的动脉的Sensitivity(表示为Sensa)高于文[20]和文[21]的方法，同时注意到，本文方法得到的对动静脉分类的4个参数结果基本一致，这是一种统计意义上的性能提高.
将由权威眼科专家确定的AVR值作为参考值，以便于与本文及其它方法进行比较. 参考值的确定包括： OD及其直径的确定、 血管类型的标注、 血管宽度的测量.
人工确定血管宽度的过程是： 在人工动静脉标定的标准图像中，选取一个长为20个像素、 宽为1个像素的“标尺”对血管的宽度进行测量. 具体地： 任意选择要测量的血管的一个位置，将标尺水平置于血管之上，数出血管内标尺的像素个数，记为N； 而后以22.5°为步长，旋转标尺，每旋转一次记录一个N值，直到标尺与图像的角度为157.5°为止，共记录8个N值； 选取这8个N值中最小的一个值作为血管的宽度.
针对何氏眼科数据库中20幅测试图像，图 7所示的Bland-Altman图分别描述了由本文方法、 另一位眼科医生以及文[10]方法所得到的AVR估计值与参考值之间的一致性. 图中，横轴为参考值与估计值之间的均值，纵轴为二者的差值，上下两条水平线代表 95％一致性的上下界，中间实线代表差值的均值. 可以看到，这3种方法的一致性界限基本一致，本文方法和另一位医生得到的AVR估计值与参考值之间的差值平均值类似，都接近于0，且几乎所有差值都比较集中地分布在一致性界限范围内，只是医生得到的差值中有一个远超过了一致性界限.
 |
| 图 7 AVR估计值与参考值之间的一致性比较Fig. 7 Agreement comparison between the estimated AVR value and AVR reference value |
表 4给出了本文方法与文[10]的AVR估计结果，以及它们与参考值之间的误差，可以看到本文方法得到的估计误差值比文[10]的结果小.
本文实现了一种眼底图像中AVR的全自动测量方法，其主要工作如下：
(1) 利用连续小波变换的多尺度特性，分别对动、 静脉图像进行了增强，提升了图像的可辨识度.
(2) 将颜色特征与结构特征相结合的特征向量作为无监督分类器的输入，提高了自动识别动、 静脉的准确度.
(3) 在感兴趣区域对动静脉进行分段，将不同子段上的AVR估计结果进行平均获得最后的估计值，减少了估计误差.
与现有一些技术方法相比，本文的方法在血管分割、 动静脉分类和AVR测量方面都具有较高的精度. AVR的测量结果与专家估计结果接近，计算时间却远少于人工方法. 因此，对基于AVR的视网膜病变早期诊断和排查系统具有一定的实用价值.
| [1] | Liew G, Sim D A, Keane P A, et al. Diabetic macular ischemia is associated with narrower retinal arterioles in patients with type 2 diabetes[J]. Acta Ophthalmological, 2015, 93(1):e45-e47. |
| [2] | McClintic B R, McClintic J I, Bisognano J D, et al. The relationship between retinal microvascular abnormalities and coronary heart disease:A review[J]. The American Journal of Medicine, 2010, 123(4):374.e1-374.e7. |
| [3] | Hubbard L D, Brothers R J, King W N, et al. Methods for evaluation of retinal microvascular abnormalities associated with hypertension/sclerosis in the atherosclerosis risk in communities study[J]. Ophthalmology, 1999(106):2269-2280. |
| [4] | Wong T Y, Knudtson M D, Klein R, et al. Computer assisted measurement of retinal vessel diameters in the beaver dam eye study:Methodology, correlation between eyes, and effect of refractive errors[J]. Ophthalmology, 2004(111):1183-1190. |
| [5] | Pakter H M, Ferlin E, Fuchs S C, et al. Measuring arteriolar-to-venous ratio in retinal photography of patients with hypertension; Development and application of a new semi-automated method[J]. American Journal of Hypertension, 2005(18):417-421. |
| [6] | Tramontan L, Grisan E, Ruggeri A. An improved system for the automatic estimation of the arteriolar-to-venular diameter ratio (AVR) in retinal images[C]// International Conference of IEEE Engineering in Medicine and Biology.Piscataway, NJ, USA:IEEE, 2008:3550-3353. |
| [7] | Tramontan L, Ruggeri A. Computer estimation of the AVR parameter in diabetic retinopathy[C]// International Federation for Medical and Biological Engineering.Piscataway, NJ, USA:IEEE, 2009:141-144. |
| [8] | Nam H S, Hwang J M, Chung H, et al. Automated measurement of retinal vessel diameters on digital fundus photographs[C]//International Federation for Medical and Biological Engineering. Piscataway, NJ, USA:IEEE, 2009:277-280. |
| [9] | Niemeijer M, Ginneken B, Abramoff M D. Automatic determination of the artery-vein ratio in retinal fundus images[C]// Proceedings of the SPIE.Bellingham, WA, USA:SPIE, 2010:7624-7628. |
| [10] | Muramatsu C, Hatanaka Y, Iwase T, et al. Automated selection of major arteries and veins for measurement of arteriolar-to-venular diameter ratio on retinal fundus images[J]. Computerized Medical Imaging and Graphics, 2011, 35(1):472-480. |
| [11] | 郭彤颖, 吴成东, 曲道奎. 小波变换理论应用进展[J]. 信息与控制, 2004, 33(1):67-71. Guo T Y, Wu C D, Qu D K. Wavelet transform theory and its application progress:A review[J]. Information and Control, 2004, 33(1):67-71. |
| [12] | Nakagawa T, Suzuki T, Hayashi Y, et al. Quantitative depth analysis of optic nerve head using stereo retinal fundus image pair[J]. Journal of Biomedical Optics, 2008, 13(6):064026-1-064026-10. |
| [13] | Muramatsu C, Hatanaka Y, Iwase T, et al. Automated detection and classification of major retinal vessels for determination of diameter ratio of arteries and veins[C]//Proceeding of SPIE Medical Imaging. Bellingham, WA, USA:SPIE, 2010, 7624:76240J1-76240J8. |
| [14] | Marin D, Manuel E G, Suero A, et al. Obtaining optic disc center and pixel region by automatic thresholding methods on morphologically processed fundus images[J]. Computer Methods and programs in Biomedicine, 2015, 2(8):173-185. |
| [15] | Otsu N. A threshold selection method from gray-scale histograms[J]. IEEE Transactions on System, Man, and Cybernetics, 1979, 9(1):62-66. |
| [16] | Ahmed M, Ward R. A rotation invariant rule-based thinning algorithm for character recognition[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2002, 24(12):1672-1678. |
| [17] | Delibasis K K, Kechriniotis A I, Tsonos C, et al. Automatic model based tracing algorithm for vessel segmentation and diameter estimation[J]. Computer Method Programs in Biomedicine, 2010, 100(2):108-122. |
| [18] | Knudtson M D, Lee K E, Hubbard L D, et al. Revised formulas for summarizing retinal vessel diameters[J]. Current Eye Research, 2003, 27(2):143-149. |
| [19] | Niemeijer M, Staal J J, Ginneken B V, et al. DRIVE:Digital retinal images for vessel extraction[DB/OL]. (2004-11-06)[2014-02-08].http://www.isi.uu.nl/Research/Database/DRIVE.. |
| [20] | Qazaleh M, Farshad T, Hamidreza P. Automated characterization of blood vessels as arteries and veins in retinal images[J]. Computerized Medical Imaging and Graphics, 2013, 37(2):607-617. |
| [21] | Dashtbozorg B, Mendonca A M. An automatic graph-based approach for artery/vein classification in retinal images[J]. IEEE Transactions on Image Processing, 2014, 23(3):1073-1083. |