第39卷  计  算  机  学  报  Vol.39 2016 论文在线出版号  No.173  CHINESE JOURNAL OF COMPUTERS  2016 论文在线出版号  No.173 ——————————————— 本课题得到国家重点基础研究发展计划项目(973 计划项目) (2014CB748600)、国家自然科学基金(81371629, 6140293, 61401294, 81401451, 81401472, 51365017)、江苏省自然科学基金(BK20140052) 、江西省自然科学基金资助项目(20132BAB203020)资助.梁礼明，男，1967年生，硕士，教授，主要研究机器视觉与医学图像分析，E-mail :lianglm67@163.com.黄朝林，男，1989年生，硕士，学生，主要研究医学图像处理的偏微分方程方法.石霏，女，1980年生，博士，讲师，主要研究方向医学图像处理及分析.吴健，男，1991年生，硕士，学生，主要研究方向机器视觉和医学图像处理.江弘九，男，1993年生，本科，学生，主要研究医学图像处理的水平集方法.陈新建(通信作者)，男，1979年生，博士，教授，主要研究模式识别、机器学习和医学图像处理与分析，E-mail : xjchen @suda.edu.cn. 融合形状先验的水平集眼底图像血管分割 梁礼明1) 黄朝林1,2)石霏2)吴健1)江弘九2)陈新建2) 1)(江西理工大学 电气工程与自动化学院，江西 赣州  341000) 2)(苏州大学 电子信息学院，江苏 苏州 215006)  摘  要  视网膜血管病变程度可作为与血管相关疾病的诊断依据，因此对视网膜血管的特征提取和分析在疾病的诊断、治疗上具有重要的临床价值。然而视网膜血管图像结构复杂，血管与背景存在灰度值交叉以及病灶和噪声的影响，致使诸多算法在相邻血管处、血管交叉处和微血管处存在分割不足的问题，且易受病灶和噪声的干扰。针对现有血管分割方法的缺陷，同时考虑视网膜图像局部血管结构与背景的可分性，以及模型的抗噪性，本文提出一种融合区域能量拟合信息和形状先验的水平集血管分割方法。在图像预处理过程中，利用GAC水平集模型控制轮廓曲线朝目标真实边界演化来获得视网膜掩模，接着通过形态学算子去除血管中心亮线以及利用高斯卷积估计视网膜背景图像，并与原图像和掩模分别进行减法和点乘运算增强视网膜血管图像；然后分析Hessian矩阵的各向异性特性，它的特征值在血管、背景和病灶上具有不同的几何性质，为最大化这种几何性质在不同结构上的差异，本文利用Hessian矩阵特征值重新构建血管响应函数，从而获得视网膜血管初步图像，作为后面的先验信息、初始化信息等，缓解水平集对初始化敏感和易误分割的问题。在水平集血管模型构建中，考虑到视网膜图像局部血管的可分性、模型的抗噪性等因素，使用RSF模型的局部区域能量拟合性质划分血管与背景，以及利用先验信息和水平集函数的几何性质构建水平集模型的形状约束项，当轮廓曲线远离视网膜血管图像的先验信息位置时就会受到惩罚，距离越远惩罚越大。完整的水平集血管分割模型包含了局部区域能量拟合项、形状约束项、面积约束项等，在能量最小化求解过程中，能够通过局部区域能量拟合性质克服现有血管分割方法在相邻血管处、血管交叉处和微血管处分割不足的问题，以及通过形状约束性质提高模型的抗噪性；最后，利用连通域面积和宽、高信息构建几何算子，进一步消除连通域面积较小的伪影和病灶，获得最终的视网膜血管分割图像。通过在以下三个数据库上的实验仿真，对HRF数据库的HRF_healthy、HRF_diabetic和HRF_glaucoma眼底图像数据，敏感度分别达到79.4547%、81.0653%、81.1773%，准确率分别达到96.1820%、94.2147%、95.6413%；对STARE、DRIVE数据库上的眼底图像数据敏感度分别达到79.0860%、75.3535%，准确率分别达到95.0340%、95.3565%。  关键词  区域能量拟合；形状先验；水平集；Hessian；响应函数；连通域   中图法分类号  TP391   论文引用格式 梁礼明，黄朝林，石霏，吴健，江弘九，陈新建，融合形状先验的水平集眼底图像血管分割，2016，Vol.39：在线出版号 No.173 LIANG Li-Ming，HUANG Chao-Lin，SHI Fei，WU Jian，JIANG Hong-Jiu，CHENXin-Jian，Retinal Vessel Segmentation Using Level Set Combined with Shape Priori，Chinese Journal of Computers,2016, Vol.39: Online Publishing No.173   网络出版时间：2016-11-25 23:11:52网络出版地址：http://www.cnki.net/kcms/detail/11.1826.TP.20161125.2311.002.html2  计  算  机  学  报                                                                                              2016年 Retinal Vessel Segmentation Using Level Set Combined with Shape Priori LIANG Li-Ming1)  HUANG Chao-Lin1,2)  SHI Fei2)  WU Jian1)  JIANG Hong-Jiu2)  CHENXin-Jian2) 1)School of Electronic and Engineering and Automation, Jiangxi University of Science and Technology, Ganzhou, Jiangxi, 341000 2)School of Electronics and Information Engineering, Soochow University, Suzhou, Jiangsu, 215000 Abstract    Thanks  to  the  severity  of  retinal  vessel  lesions  is  an  important  reference  to  the  diagnosis  of  vessel  diseases;  the segmentation and analysis of the retinal vessel have a great clinic value in diagnosing and treating these diseases. However, due to complex structure of the retinal blood vessel, gray level crossing and the influences of the focus and noise, many algorithms have disadvantage of inadequate segmentation on the position of adjacent vessels, crossing vessels and micrangium as well as being affected easily by the focus and noise. Aiming to the shortcomings of existed methods and considering that partial vessel structure can be separated from background and the resistance of model to noise, we propose level set method to segment blood vessel fusing regional energy fitting information and shape priori. First, in the stage of preprocessing, it uses GAC level set model to control contour curve evolving to real border, which leads to retinal mask, eliminates central bright line of vessel and evaluates background image of retina, conducts subtraction and multiplication with initial image and mask to reinforce retinal vessel image. Than by the characteristic analysis of anisotropy of Hessian, its character values have different geometric property on vessels, background and the focus, to maximize the difference on different structures, we rebuild vessel response function by the eigvalues of  Hessian  matrix  to  acquire  primary  image  of  retina  as  prior and  initialization  information  thereafter,  which  alleviate  the sensitivity  of  level  set  model  to  initialization  and  missegmentation.  On  constructing  the  model  of  level  set,  taking  into consideration separation of local retinal image, noise-resistance of model and so forth, we make a distinction between vessels and background  via  the  property  of  regional  energy  fitting  of  RSF  model  and  build  shape  constraints  of  level  set  using  prior information  as  well  as  geometric  property  of  level  set  function,  the  curve will  receive  punishment  when  it  is  away  from  the position of prior information, the further, the bigger the punishment is. The integral vessel segmentation model of level set is composed of local energy fitting, shape constraint, area constraint; on the process of minimization, it overcome the shortcomings by  energy  fitting  of  insufficient  segmentation  at  adjacent  vessels, microvascular  and  vascular  intersection  that  exist  in  many established method, and improve the noise-resistance of model via shape constraint. At last, the final retinal segmentation image is achieved by area of connected domain and geometric operator constructed by width and height information which can further eliminate the fake shadow of connected domain and the focus. By the experiments on following three databases, the sensitivity is 79.4547%、81.0653%、81.1773%, accuracy is 96.1820%、94.2147%、95.6413% corresponding to dataset of fundus image of HRF_healthy、HRF_diabetic  and  HRF_glaucoma  in  database  of  HRF;  on  STARE、DRIVE, the  sensitivity is 79.0860%、75.3535%,accuracy  is  95.0340%、95.3565%  respectively.  The  results  above  show  the  effectiveness  of  vessel  segmentation method proposed in this paper.  Key words    regional energy fitting ; shape priori; level set; Hessian; response functions; connected domain  1 引言 中国是世界上盲人数量最多的国家，约占全球盲人总数的18%[1]，其中视网膜血管病变是致盲的重要原因之一。青光眼、老年性黄斑、高血压、动脉硬化、糖尿病、心血管疾病等全身性疾病[2,3]均可引起视网膜血管病变。而在眼科学领域，视网膜血管数量、分支、角度、宽度等都可作为与视网膜血管相关疾病的诊断依据[4]。因此，通过数字图像处理技术分割视网膜血管，做定性和定量分析，可辅助医师诊断、分析患者病情。做到早发现、早治疗，从而降低视网膜疾病患者的视网膜损伤和致盲概率。2D眼底图像可利用眼底相机直接拍摄得到，且能清晰地反映视网膜上的微血管和其它病灶[5]，因此具有重大的医学应用价值。 目前在国内外的血管分割算法有多尺度血管增强滤波[6-9]、多阈值的血管检测[10]、以形态学为基础技术血管分割[11,12]、线性检测[13,14]、使用神经网络的血管分割算论文在线出版号  No.173                          梁礼明等：融合形状先验的水平集眼底图像血管分割                                      3  法[15-18]、多尺度层分解和局部自适应阈值血管分割方法[19]以及基于活动轮廓模型的血管分割[20-23]等。其中文献[6-9]通过分析Hessian矩阵在管状结构和非管状结构上的特征值差异性，或进一步利用特征值构建血管响应函数估计血管图像，但因为高斯卷积算子的影响，分割结果存在相邻血管易相连、血管过宽、微血管和血管交叉处分割不足等问题，导致用于复杂血管结构的视网膜血管分割时精度较低。文献[10-12]以阈值和形态学等技术为基础的血管分割方法，存在对图像噪声过于敏感和难以克服目标与背景存在灰度值交叉等问题。文献[13,14]分析视网膜血管图像局部性质，沿着血管方向上的检测线的像素均值大于其它方向上检测线的像素均值，以此增强血管图像，但在血交叉处易存在分割断裂和对病灶产生误分割等问题。文献[15-18]使用大量的血管图像训练神经网络，把血管分割作为二分类问题，对噪声具有较强的适应性，但依然无法克服血管与背景灰度值交叉和微血管易断裂问题。文献[19]使用多尺度小波核匹配滤波器增强视网膜血管图像，然后使用多尺度层分解和局部自适应阈值获得血管的二值图谱，解决了相邻血管分割时易相连的问题，但对噪声较为敏感。文献[20-22]提出基于活动轮廓模型的视网膜血管分割算法，利用轮廓线可在曲线自生内力和图像数据的外力驱动下自动实现拓扑结构变化，从而在轮廓线的进化过程中逼近真实的血管边缘，克服了相邻血管相连、血管过宽、血管交叉处分割不足等问题，但都不同程度地存在对视网膜微血管分割不足和对噪声、病灶等敏感等问题。文献[23]虽然利用空间模糊聚类性质建立了一个新的活动轮廓模型，成功解决对微血管、血管交叉处分割不足等问题，但对初始化较为敏感，导致无法直接用于具有复杂拓扑结构的视网膜血管图像分割，以及时间开销大。 活动轮廓模型是医学影像图像分割的经典方法之一，分为参数主动轮廓模型[24-26]和几何主动轮廓模型[27-36]。其中几何主动轮廓模型的出现源于S.  Osher等人在文献[37]中成功把水平集函数嵌入能量泛函中，因此又俗称水平集，它把N维曲线或曲面的演化隐藏于N+1维曲面或立体空间中，可在无先验信息下实现拓扑结构的任意变化。水平集的特点无疑非常适应于具有复杂结构的医学影像图像分割，因此在医学图像分割领域中得到飞速发展。但不论经典的测地线活动轮廓(Geodesic active  contours，GAC)模型[30]，还是基于区域特征的C-V(Chan T F and Vese L)模型[36]都有致命缺点，在迭代过程中都需周期性地使用符号函数重新初始化水平集函数，维持它的平滑性，具有时间开销大以及可能造成分割失败等问题。对此，Li在文献[31]中通过构建关于扩散模型的积分表达式，并融入到能量泛函中，使水平集函数在迭代过程中依靠内部能量扩散项逐渐逼近符号距离函数，从而自动实现水平集函数的光滑性维持。但依靠图像梯度信息的GAC模型在消除曲率为负的凹陷处分割时易造成血管弱边缘泄漏。C-V模型成功解决GAC模型易造成弱边缘泄漏和凹陷处难分割的问题，但它把目标与背景简单的考虑为全域二聚类问题，因此难以实现血管与背景存在灰度值交叉的视网膜血管图像分割。Li在文献[32]中引入高斯模板函数，把C-V模型的二聚类转化为曲线内、外部两个局部区域能量拟合(region-scalable fitting energy，RSF)问题，背景与目标存在灰度值交叉的分割问题得以解决。但又因为水平集对初始化敏感和噪声敏感，因此RSF模型依然无法直接用于视网膜血管图像的分割。 针对现有方法在视网膜血管分割上的不足，本文首先利用形态学算子和高斯卷积算子增强视网膜图像；然后利用Hessian各向异性特性和改进的血管响应函数粗略地分割视网膜血管图像，并以此作为形状先验、初始化信息等，提出一个融合区域能量拟合信息、形状先验的水平集模型，从而成功分割视网膜血管图像；最后利用连通域面积Area和宽、高( ) , WH构建几何算子，消除连通域面积较小的伪影和病灶，获得最终的视网膜血管分割图像。 2 眼底图像预处理 本文使用HRF、STARE、DRIVE数据库的眼底图像作为实验数据。通过对RGB眼底图像绿色、红色、蓝色通道图像对比，发现绿色通道具有较多的血管信息和较高的对比度，因此取眼底图像的绿色通道图像作为血管分割的实验数据。为使分割更为精确，采用文献[30]中的GAC模型自动获取每个眼底视网膜图像“掩模”，并对眼底视网膜图像做边缘扩展；然后估计眼底视网膜图像背景，并与原图像做减法运算和灰度值拉伸，获得灰度值在一个固定范围视网膜血管增强图像。       2.1 获取视网膜掩模 取绿色通道图像I，然后利用文献[30]的GAC模型获得视网膜“掩模”。因为曲线可在内力和基于梯度信息的外力作用下自动向视网膜真实边界逼近，因此可以自动获得高精度的视网膜“掩模”，这将有利于做视网膜边界扩展，且在后续使用卷积与减法运算增强血管图像时，能够获得更加光滑的视网膜图像边界。 图1展示了HRF数据库正常人的第12_h个眼底图像视网膜“掩模”的获取过程，图像尺寸从2336×3504降4  计  算  机  学  报                                                                                              2016年 维到1176×1752像素。(a) 、(b) 、(c)、(d)分别表示绿色通道图像、分割初始化、分割结果、“掩模”。  (e)、(f)分别表示(b)、(e)图像中白色小方框部分的细节对比。  (a)绿色通道图像       (b)初始化                    (c)分割结果           (d)掩模      (e)来至(b)的局部放大 (f)来至(c)的局部放大 图1  使用GAC模型自动获取视网膜“掩模” 2.2 血管图像增强     视网膜图像因为噪声和光照变化使血管图像对比度较低，影响血管分割算法性能。这里对视网膜图像作边缘扩展、血管中心亮线消除、背景估计、灰度拉伸等处理，从而增强视网膜血管图像。对HRF数据集，视网膜图像统一从2336×3504降维到1176×1752像素，其中血管中心亮线直径在3个像素左右，因此背景估计时的高斯模板尺寸为100×100，血管中心亮线消除时形态学算子模板尺寸为5×5；这里对STARE、DRIVE数据集，高斯模板尺寸为50×50，形态学算子模板尺寸为3×3。血管增强具体步骤如下： 第一：利用上面获得的视网膜“掩模”信息，对图1(a)做基于镜像对称的边缘扩展，边缘扩展的尺寸等于“第二”中的高斯模板尺寸。 接着使用板尺寸为5×5形态学算子，对图像作开运算，消除血管中心亮线，防止含有区域能量拟合项的水平集模型产生误分割。ExtenI是消除了血管中心亮线的视网膜边缘扩展图像，如图2(a)。 第二：构建尺寸为100×100的高斯模板，对图像ExtenI做卷积运算，估计视网膜图像背景BackgI，并分别与原图像、“掩模”做减法和点乘运算，从而获得增强的视网膜血管图像enhanceI。 第三：把图像enhanceI的灰度值线性拉伸到[0-255]，获得图像GraySI。对血管增强后图像做灰度拉伸是因为包含区域能量拟合项的水平模型，在能量函数最小化迭代求解过程中，能量并不趋于零，因此目标和背景灰度值差需在一个较大的动态水平值上。 本文血管增强方法不同于文献[15]。本文在对视网膜做边缘扩展基础上采用高斯卷积估计视网膜背景，获得了一个背景更加均匀的血管增强图像，是因为边缘扩展能够消除视网膜原始边界处的灰度突变。下面图2展示了使用本文方法对第12_h眼底图像的视网膜血管增强过程，(a)、(b)、(c)、(d) 分别表示视网膜边缘扩展、视网膜背景估计、血管增强、灰度值拉伸的结果。从图(c)和(d)可以看出，视网膜血管图像得到增强，黄斑中心和视盘中心灰度值的剧烈变化得到抑制。      (a)ExtenI边缘扩展        (b)BackgI背景估计      (c)enhanceI血管增强       (d)GraySI灰度值拉伸 图2  视网膜血管图像增强 3 视网膜血管分割方法 本文在分析Hessian矩阵各向异性及特征值、特征向量基础上，改进血管响应函数，利用血管响应函数和阈值得到血管的粗分割图像，并以此信息初始化水平集函数和构建形状先验项等，它克服了水平集模型对初始化敏感以及区域能量拟合项对噪声敏感的问题。此外，利用连通域面积Area和宽、高( ) , WH构建几何算子，消除水平集模型分割结果中的病灶和连通域面积较小的伪影，获得最终的视网膜血管分割图像。 3.1 初步估计视网膜血管 3.1.1 Hessian矩阵的几何意义分析 一阶导数和二阶导数能够描述物体局部的特征变化，因此在基于偏微分方程的数值图像处理上常用于构建能量驱动项和特征描述项等。例如，可用含一、二阶导数的梯度构建水平集模型边缘项的方向指示函数，可用二阶导数描述线状、管状等物体特征。 考虑到眼底视网膜血管图像恰好呈现出一种管状结构，对此本文引入由二阶偏导数构建的Hessian矩阵。又因为数值图像常存在噪声和伪影，一般将原图像与尺度为s的高斯函数的二阶偏导数做卷积运算得到Hessian矩阵。Hessian矩阵的定义如下： 论文在线出版号  No.173                          梁礼明等：融合形状先验的水平集眼底图像血管分割                                      5  ( )( ) ( )( ) ( )222222,,,*,, GraySG x y G x yy x yH x y IG x y G x yy x ysssséù¶¶ êú¶ ¶ ¶êú =êú¶¶ êú¶ ¶ ¶ êú ëû     (1) 二维Hessian矩阵为实对称矩阵，有两个特征值1l和2l。假设在视网膜血管图像中，1l对应的特征向量表示血管的走向，2l对应的特征向量表示与血管走向垂直的方向，那么12<< ll，即在血管上的一点( ) , xy，其走向方向的曲率远远小于走向垂直方向的曲率，如图3。而在非血管的背景区域，因为无明显的管状结构，各点上的特征值有12 ll»，若12>>0 ll»，那么表示此点与周围点存在明显的灰度值变化(为噪声点或者病灶)；若120 ll»»，那么表示此点与周围点不存在明显的灰度值变化。20 l <表示血管比背景亮(为凸)，20 l >表示血管比背景暗(为凹)。本文的视网膜血管图像，应取20 l >。  图3  血管上的特征值描述 3.1.2 构建血管响应函数         根据Hessian矩阵对图像不同点( ) , xy的特征值描述，就可构建用以区别血管与背景的描述算子。Frangi等[8]提出使用“脊选评分”  (Ridgness score)21 1 / / Rbll =描述血管，以及用Hessian矩阵的F范数去除背景和微弱信号的影响。二阶Hessian矩阵的F范数定义如下： 221iFiSH l===å                      (2) 因为背景中非噪声点灰度值在相互垂直的两个方向都不存在突变， 0, 1, 2ii l »=，F范数和Rb较小。而在血管上因12<< ll和20 l >，F范数和Rb都较大，因此可成功划分血管与背景。但在噪声点12>>0 ll»，Rb虽依然较小，但F范数较大，因此难免产生错误的血管响应。对此本文构建了一种新的血管描述算子，这里定义为J。 ( ) ( )22 1 2 1 2- J l l l l =+                (3) 因此，当点( ) , xy在血管上时，因12<< ll，J较大。 当点( ) , xy在噪声或背景上时，因12 ll»，J较小。 根据以上性质，本文构造血管响应函数( ) , xysn用以区分血管与背景，从而初步估计眼底视网膜血管图像。( ) , xysn定义如下： ( )222 2max00 ,1 exp 1 exp 02R xy JJb slnlb< ìïæö æöæö æö =íç÷ - - - - >= ç÷ç÷ ç÷ ïç÷ ç÷ ç÷ç÷ èø èø ïèøèø î  (4) 其中b是尺寸参数，通常取0.5 b=，s是高斯核函数尺寸，maxJ是J的最大值。这样式(4)能最大化血管像素的响应，使血管响应函数 ( ) ,1 xysn ®，而在非血管背景平坦区域上使血管响应函数( ) ,0 xysn ®，即( ) [ ] , 0,1 xysn Î，这是本文血管相应函数与文献[8]的区别，如图4(a)和(b)。因为眼底视网膜血管图像尺寸不一，因此需在多尺度下计算血管响应函数( ) , xysn，最后取各尺度下的最大响应值，其定义如下： ( ) ( )min max, max , x y x yss s snn££=                     (5) 其中mins和maxs为感兴趣的眼底视网膜血管的最小和最大尺寸。对HRF数据集[ ] 1,11 sÎ，对HRF数据集STRAE和DRIVE数据集[ ] 1, 7 sÎ，步长为2。通过式(4)和(5)就能够在增强后的图像GraySI上计算( ) , xy n，获得血管响应图像RespI。 根据( ) [ ] , 0,1 xysn Î，对图像RespI进行二值化。取较小的阈值10.00001 c =获得二值化图像1BinaryI，作为水平集函数j的初始化，如图4(c)。取略大的阈值20.0001 c =获得二值化图像2BinaryI，并使用尺寸为2×2的模板对2BinaryI作开运算和删除连通域面积小于100 的伪影(参考文献[14])，得二值图像B I，如图4(d)。这样能使初始化中包含更多的血管信息，在使用B I构建形状约束时含有较少的非血管结构，从而提高模型的抗噪能力和血管分割精度。     (a)文献[8]的血管响应           (b)RespI 血管响应     (c) 1BinaryI                    (d)B I二值图像 6  计  算  机  学  报                                                                                              2016年 图4   血管响应函数响应结果 3.2 建立视网膜血管分割的水平集模型 使用视网膜血管粗分割结果，就可构建一个包含区域能量拟合项、形状约束项、正则项、长度惩罚项、面积约束项的水平集模型。其中，区域能量拟合项具有局部区域数据拟合能力，它可克服背景与血管存在灰度值交叉的问题，保证相邻血管、血管交叉处和微血管处的正确分割；形状约束项和面积约束项可克服区域能量拟合项噪声敏感的问题；正则惩罚项可使水平集函数在演变过程中具有正则性，保证水平集函数在演变过程中获得精确的数字解；长度惩罚项可使分割轮廓曲线具有光滑性。本文利用区域能量拟合项、形状约束项、正则项、长度惩罚项、面积约束项的相互作用，完成在血管相邻处、血管交叉处和微血管处的正确分割，同时克服噪声干扰。其能量模型定义如下： ( ) ( ) ( )( ) ( ) ( )1 2 1 2, , , ,RSR P RE f f E f f EP L Af f h fm f n f g f=++ + +     (6) 其中式(6)的右边分别为局部区域能量拟合项、形状约束项、水平集函数正则性维持项、长度惩罚项、面积约束项。, , , h m n g为权重系数。  3.2.1 区域能量拟合项 RSF两项水平集模型利用曲线内、外部两个函数( ), 1, 2if x i =进行局部区域灰度值拟合，其局部能力主要通过高斯核函数( ) k x ys-加以限制。以y为中心，当x逐渐远离y时，表示x点能量拟合权重越小，因此可克服视网膜血管与背景存在灰度值交叉的问题，从而引导分割曲线向血管边缘逼近。 区域能量拟合项定义如下： ( )( ) ( ) ( ) ( ) ( ) ( )12221,, Ri i iiE f fk x y I y f x M y dy dxesfbf ==-- åòò  (7) 其中s为高斯核尺寸参数，, 1, 2ii b =为区域能量拟合项权重系数，f为Lipschitz水平集函数，( ), 1, 2iMief =为区域函数，表示曲线内部和外部。定义( )1() MHeeff=，( )2( ) 1 MHeeff=-。( ) Hef为具有正则性的Heaviside函数，其导数为Dirac函数( )edf。 ( )121 arctan2He ffpeæö =+ç÷èø èø           (8) ( )221eedfp e f=+                         (9) 其中e为尺寸参数，一般取1 e=。         针对视网膜血管分割，本文取3.1.2节中1BinaryI作为水 平 集 函 数 的 初 始 化 信 息 ， 即( ) ( )1Binary, , 0 2 1 2 x y I f =-。 3.2.2 形状约束项 基于区域能量拟合项( )12 ,, RE f f f的视网膜血管分割对噪声敏感，易把部分非血管的噪声区域错分割为血管，造成精度丢失。为克服区域能量拟合项的不足，本文利用3.1.2 节获得二值图像B I构建视网膜血管的形状约束项。 首先根据二值血管图像B I初始化如下水平集函数： ( )20 , , 020BBIxyIf= ì=í-¹ î                 (10) 然后使用GAC模型[30]或RSF模型[32]控制轮廓曲线C的演变，在时刻*t获得符号距离函数( )*, xt j表示视网膜血管形状： ( )( )( )*,,0,outsideinsided x C x Cx t x Cd x C x CjÎ ìï=Îíï-Î î      (11) 其中，( ) , d x C表示视网膜图像中点x到轮廓曲线C的Euclidean距离，outsidexC Î表示点在轮廓线C外面，insidexC Î表示点在轮廓线C内部，xC Î表示点在轮廓线C上。         在水平集分割模型中，通过最小化能量表达式控制轮廓曲线C朝目标边缘演化。因此，可令( )*, xt j表示视网膜血管先验形状，构建关于水平集函数f和( )*, xt j的能量泛函表达式( )SE f。 形状约束项能量表达式( )SE j定义为： ( ) ( ) ( )2SE H H dxee f f jW=-ò           (12) 在最小化能量泛函过程中，当视网膜血管分割轮廓曲线偏离形状约束先验时就会受到惩罚，偏离距离越大惩罚越大，从而克服区域能量拟合项对噪声的敏感。 3.2.3 其它能量项 水平集函数f在进化过程中必须始终具有正则性，这种正则性要求水平集函数不能太平坦或者太陡(至少在零水平集附近)，才能保证水平集函数平稳地进化，从而控制轮廓曲线C逐渐逼近目标的边缘，获得精确的数值解。因此本文引入文献[21]中的水平集正则化项( )RP f，其定义如下： ( ) ( )2 112RP dx ffW= Ñ -ò                 (13) 在能量最小化求解过程中，为使轮廓曲线C尽可能的短和光滑，以及消除由视网膜图像中的噪声和病灶引起的局部最小值问题，还需在能量表达式中，添加轮廓论文在线出版号  No.173                          梁礼明等：融合形状先验的水平集眼底图像血管分割                                      7  曲线C的长度惩罚项( )PL f，其定义如下： ( ) ( ) PL H dxeffW=Ñò              (14) 为进一步克服视网膜血管图像噪声和病灶引起的误分割，本文利用二值图像B I点乘预处理图像GraySI，获得血管图像ReaI。然后对ReaI高斯卷积求一阶偏导，获得梯度图像ReaI Ñ。以此构建融合血管边缘梯度信息的加权面积约束项( )RA f。图像ReaI和ReaI Ñ如图5(a)和(b)所示。 加权面积约束项( )RA f的定义如下： ( ) ( ) ( ) 1RA g H dxeffW=-ò                (15) 其中( )2Re1 / 1 *ag G I = + Ñ为非负单调递减函数，G为高斯函数。根据式(15)的定义，在非血管的背景区域1 g=，面积约束项具有一个较大的惩罚，而在血管附近0 g®，面积约束项只获得一个较小的惩罚。     (a)ReaI                               (b)ReaI Ñ 图5  血管和血管的梯度图像 3.3 移除病灶和伪影 在HRF、STARE、DRIVE数据库上，使用本文的水平集模型分割出的视网膜血管图像任然含有部分连通域面积较小的伪影和病灶。这里构建8连通域模板标记水平集分割结果，从而获得每个连通域的面积Area和宽、高( ) , WH信息。如果一个连通域的像素是非血管像素，那么它的连通域的宽高应相近，且连通域面积和宽高乘积大小相近，反之，若是血管像素，那么是狭长状。因此/ WH和WH´必须同时满足以下两个性质： 1、0.4 / 2.5 WH << 2、3.5 W H Area ´< 这里同时定义8连通域面积小于30为伪影： 3、30 Area<。 利用几何算子1和2、3能够移除伪影和部分病灶，提高视网膜血管的分割精度，其结果将在实验5.3给予展示。 4 能量函数求解和参数设置 4.1能量函数求解 本文使用标准梯度下降流法求解能量函数。首先固定水平集函数f，最小化关于变量( ), 1, 2if x i =的能量函数( )12 ,, E f f f。因此获得拟合函数值( ), 1, 2if x i =可通过计算满足如下关系的欧拉方程： ( ) ( ) ( ) ( ) ( ) ( ) 0ii k x x I y f x M y dyes f - - =ò  (16) 求解式(16)可以获得： ( )( ) ( ) ( ) ( )( ) ( ) ( )*, 1, 2*iiik x M x I xf x ik x M xesesfféù ëû ==  (17) 然后固定( ), 1, 2if x i =，最小化关于变量f能量函数( )12, E f f f。能量函数最小化求解可用如下水平集函数的梯度下降 ( )( )( ) ( )( ) ( ) ( ) ( )( )1 1 2 221eetg divHHdiveeee e efd f b bfg d f nd ffhd f f jfmff¶= - -¶æöÑ++ ç÷Ñèø--æöÑ+ Ñ - ç÷Ñèø      (18) 其中( )edf是来至式(9)的Dirac  delta函数。ie是满足如下性质表达式： ( ) ( ) ( ) ( )2, 1, 2ii e x k x y I x f y dy is = - - =ò   (19) 因水平集正则项( )RP f在能量最小化过程中有自动驱使1 f Ñ®的性质，所以可对式(18)左、右边分别采用前向差分和中心差分的数值离散化方案求解。 前向差分：           ( )1:nnn niitiiDttff ff++ - ¶»= ¶D            (20) 中心差分：          ( )( )011110:2:2nnn niixiinnn njjyjjDxxDyyff ffff ff+-+-ì- ¶»= ï¶D ïí-¶ï»= ï¶D î        (21) 4.2 参数设置 通过血管增强和灰度值拉伸等数值图像处理方法获得的视网膜图像GraySI依然存在噪声和病灶，在各能量项权重系数都设为固定值情况下，不同的噪声水平和病灶量决定血管分割有不同的精度和时间开销，这就意味着固定的权重系数设置不能适应分割所有的视网膜血管图像。虽然可采用高斯卷积等滤波方案进一步平滑图像GraySI，但会同时模糊视网膜的微血管，造成对微血管的分割不足，引起分割精度丢失。 为获得满意的血管分割结果和固定时间开销，本文8  计  算  机  学  报                                                                                              2016年 首先利用视网膜血管图像GraySI与3.2.3节中ReaI图像做减法运算，获得粗略去除血管信息的视网膜背景图像RvbI，如图6(a)，(a)中的黑色线条表示血管信息已被粗略移除。计算RvbI的方差2s，利用方差信息设置长度惩罚项、面积约束项、形状约束项的权重系数，从而固定水平集模型分割视网膜血管的迭代次数为40次，以及提高本文算法在不同数据库上的适应性。  (a)RvbI  图6  不含血管的背景图像 长度惩罚项系数： ms n=                           (22) 面积约束项系数： logans g =                    (23) 形状约束项系数：  2ls h=                           (24) 式(22)和式(23)惩罚系数不易过大，特别是面积约束像系数，它必须是一个缓慢变化的值，否则会惩罚过度，引起分割精度丢失；形状约束项用于克服数据项对噪声和病灶敏感的问题，因此需要一个较大的值。这样在血管分割时既可克服噪声和病灶等因素干扰，又能固定时间开销。根据大量实验，本文m n a l 、、、分别为5、1、7、1。 5 实验结果与分析 实验时使用的仿真平台为Matlab  R2012b,计算机配置为Intel(R) Core(TM)i7-3770 CPU @ 3.40GHz,8.00G内存,64-bit Operating System。本文采用中心差分和前向差分对水平集模型进行离散化数字求解，所以必须满足1 / 4 tm<[31]，这样才能保证水平集演变稳定，因此在所有的血管分割实验中，时间步长0.1 t =，正则化项系数1 m=。区域能量拟合项系数121 bb==，长度惩罚项、面积约束项、形状约束项分别按照式(22)、(23)、(24)设置系数,, n g h。分割一个视网膜血管图像的时耗：在HRF数据库上是1分钟左右，在STARE 和DRIVE上约20秒。 为证明本文算法的血管分割性能，分别在HRF、STARE、DRIVE数据库上进行大量的血管分割实验仿真。其中，HRF(high-resolution fundus image database)于2013年由捷克Brno科技大学和Erlangen-Nurembeg大学合作公开发布，分别包括15健康人、糖尿病和青光眼患者眼底图像，图像尺寸2336×3504，在血管分割时统一降维到1176×1752像素；STARE(structured snalysis of the retina)于1975年由Michael  Goldbaum建立，在2000年公开，它分别包括10幅有病变的眼底图像和10没病变的眼底图像，图像尺寸605×700；DRIVE(digital  retinal images for vessel extraction)于2004年由Niemeijer团队建立，分别包括7幅有糖尿病视网膜病变的眼底图像和33幅没糖尿病视网膜病变的眼底图像，它给出了20幅手动血管分割金标准结果，图像尺寸768×584。此外，在血管分割实验中取包含较多血管信息和较高的对比度的绿色通道图像作为血管分割图像。 5.1 视网膜血管增强算法对比 本文对视网膜图像进行了边缘扩展，且采用模板为100×100的高斯卷积模糊视网膜血管图像，它等价于使用线性扩散模型。下面图7(a)和(b)展示了文献[15]和本文的视网膜血管增强结果，可从白色线之间的区域及局部放大图(c)和(d)看出，本文算法获得的视网膜增强图像边界灰度值更加均匀和平滑。此外，选用在4.2 节中粗略去除血管信息的视网膜背景图像RvbI的方差2s作为度量背景平滑性和噪声水平的评价因子。实验数据使用HRF数据库的眼底图像，结果展示在下面表1中。HRF_healthy和HRF_diabetic分别表示健康人和糖尿病患者，“h”和“dr”表示不同的眼底图像。     (a)文献[15]                                            (b)本文方法      (c) (a)的局部放大                (d) (b)的局部放大      图7  血管增强结果对比 从下面表1可知，无论是健康人还是糖尿病患者的眼底图像，本文算法获得的血管增强图像，其背景灰度值间具有更小的方差，证明了本文获得的视网膜血管增强图像噪声水平更低和更具有平滑性。 表1  本文算法与文献[15]获得的血管增强图像方差对比 HRF_healthy  01_h  03_h  05_h  07_h  09_h 本文方法2s  71.9248  95.08575  87.6188  57.3860  59.7692 论文在线出版号  No.173                          梁礼明等：融合形状先验的水平集眼底图像血管分割                                      9  文献[15]2s  106.3369  112.8108  157.8549  91.5547  101.1754 HRF_diabetic  01_dr  03_dr  05_dr  07_dr  09_dr 本文方法2s  73.4818  175.9325  70.3389  82.4704  139.0179 文献[15]2s  177.9209  284.7996  77.3721  168.5368  191.3803 5.2 视网膜血管分割结果 这里分别展示来至HRF数据库的健康人(HRF_healthy)、糖尿病患者(HRF_diabetic)和青光眼患者(HRF_glaucoma)的第12_h 、01_dr、06_g个眼底图像的分割结果，如图8，图像尺寸统一从2336×3504降维到1176×1752像素。  (a) 12_h绿色通道图像  (b) 12_h自动分割   (c) 12_h手动分割  (d) 01_dr绿色通道图像 (e) 01_dr自动分割   (f) 01_dr手动分割  (g) 06_g绿色通道图像  (h) 06_g自动分割    (i) 06_g手动分割 图8  HRF数据库的视网膜血管分割结果 图8展示了来至HRF数据的眼底视网膜血管图像的分割结果。其中(a)、(d)、(g)分别表示健康人、糖尿病患者和青光眼患者的绿色通道图像，(b)、(e)、(h)表示使用本文算法与之对应的自动分割结果，(c)、(f)、(i)是与之对应的某位专家的手动分割结果(金标准)。从(a)、(b)、(c)可知本文算法能最大限度地提取健康人的视网膜血管图像，分割结果接近金标准。从(d)、(e)、(f)可知，尽管糖尿病患者视网膜中有大量的病灶，但本文算法依然能较大程度的克服病灶的干扰，降低误分割率，获得一个较好的分割结果。从(g)、(h)、(i)可知，尽管青光眼患者视网膜模糊，造成血管不够清晰(特别是微血管)，但从分割结果(h)可知，依然获得一个较好的分割结果。以上分割结果的定性展示证明了本文算法在正常和病变视网膜血管图像分割上具有较强的鲁棒性和较高的精度。 5.3 移除伪影和病灶 自动分割结果存在连通域面积较小的伪影和病灶，假阳性值(PF)较高，影响分割精度。这里利用3.3中的3个几何算子去移除伪影与病灶。图9展示了移除图8(b)、(e)、(h)中的伪影和病灶，获得最后的血管图像分割结果。  (a)12_h伪影和病灶  (b)01_dr 伪影和病灶 (c)06_g伪影和病灶  (d) 12_h分割结果    (e)01_dr分割结果    (f) 06_g分割结果 图9   移除伪影和病灶 从图9(a)、(b)、(c)可知，连通域面积较小的伪影成功地被移除，获得了最终的分割结果(d)、(e)、(f)。特别是糖尿病患者(HRF_diabetic)的眼底图像，因为含有大量的病灶，即使包含区域能量拟合项和形状约束项的水平集血管分割模型，也难以一次性的排除病灶干扰。而从(b)可知，通过3个几何算子，大量病灶得以移除。这里证明了这3个几何算子对移除连通域面积较小的伪影和病灶的有效性。 5.4 分割算法对比 为证明本文融合形状先验的水平集血管分割模型的优越性，分别与多尺度血管增强滤波[8]、线性检测[14]、GAC模型[30]和C-V模型[36]做对比(这里使用Hessian矩阵血管粗分割结果作为GAC和C-V模型水平集函数的初始化)，最后展示了某位专家手动分割的金标准。实验数据使用HRF数据库健康人的第12_h个2D眼底图像，维数从2336×3504降维到1176×1752像素，且把不同算法得到的二值血管图像点乘去除中心亮线的绿色通道图像，以此作为分割结果的比较，这样可清晰、直观地看出分割出的视网膜血管尺寸是否偏大。图10中，第二、三列分别对应左边第一列中两个白色方框部分，“暗色”的管状结构表示血管，“白色边缘”表示误分割，白色字母“A”、“B”和“C”分别表示相邻血管、血管交叉位置和微血管。 10  计  算  机  学  报                                                                                              2016年        (a)文献[8]Hessian                                     (b)                                                     (c)        (d)文献[14]线性检测                       (e)                                                   (f)          (g)文献[30]GAC                                          (h)                                               (i)                                            (j)文献[36]C-V                                      (k)                                                 (l)        (m)本文算法                                               (n)                                                   (o)        (p)手动分割(金标准)                                     (q)                                               (r)                                                                  图10   比较不同算法的视网膜血管分割结果                   论文在线出版号  No.173                          梁礼明等：融合形状先验的水平集眼底图像血管分割                                    11  从图10第一列(d)、(g)、(j)与(m)比较可知，本文方法比文献[14]、文献[30]、文献[36]、分割出的视网膜血管数量更多，与文献[8]分割出的血管数量相近。但从(b)和(c)知，文献[8]基于Hessian矩阵和血管响应函数分割出的视网膜血管尺寸明显大于血管真实尺寸，具有较高的假阳性值。白色字母“A”、“B”和“C”展示了不同算法在相邻血管处、血管交叉处和微血管处的分割性能。其中，从“A”知除本文算法和文献[14]之外，文献[8]、文献[30]、文献[36]均不能正确地完成相邻血管的分割；从图10 (e) 中“B”知文献[14]因血管交叉处的血管尺寸接近最大模板尺寸，致使分割出的血管产生断裂，文献[8]、文献[30]、文献[36]在血管交叉处的凹陷位置把大量的非血管像素误分割为血管像素，但本文算法克服了在“B”处其它文献的不足；从“C”知对于微血管结构，文献[14]和文献[36]无法完全检测出其血管结构，其上图10 (e)和(k)展示了文献[14]和文献[36]在“C”处存在严重的血管漏检，而从10图(n)、(o)可看出本文算法能够最大限度地检测出微血管结构，且成功地克服了部分微血管易断裂问题。与金标准(p)、(q)和(r)比较可知，本文算法因同时克服了在相邻血管处、血管交叉处和微血管处其它算法的不足，使分割出的血管结构最接近金标准和血管真实尺寸值。 5.5 不同算法分割结果的定量分析 为进一步证明本文算法的视网膜血管分割性能，下面将在HRF和STARE、DRIVE数据库上进行大规模的视网膜血管分割实验，并作定量分析。让本文方法分别与文献[8]、文献[9]、文献[14]、文献[15]、文献[19]、文献[20]、文献[22]、文献[30]和文献[36]算法做对比，其中文献[8]、文献[14]、文献[30]和文献[36]是作者本人编程实现，其它对比文献是直接截取其论文结果。 这里定义如下几个视网膜血管分割性能评价指标： 敏感度（Sensitivity，又称真阳性率）：正确分割血管像素的总数与真实血管像素的总数之比。 假阳性率（FPF）：错分背景像素为血管像素的总数与真实背景像素的总数之比。 特异性（Specificity）：正确分割背景像素的总数与真实背景像素的总数之比。 准确率（Accuracy）：正确分割血管和背景像素的总数之和与整个图像总像素之比。 F-measure：敏感度与精确率(Precision)之间的比率关系。 相应的计算公式定义如下： PPNTSensitivityTF=+                    (25) PPNFFPFFT=+                         (26) NNPTSpecifityTF=+                         (27) rPPPTP ecisionTF=+                       (28) PNP N P NTTAccuracyT T F F+=+ + +                  (29) 2 PrecisionSensitivitySensitivity PrecisionF measure - = ´´+  (30)其中, , ,P N P NT T F F分别表示真阳性、真阴性、假阳性、假阴性。真阳性（PT）：正确识别血管像素为血管像素的数目。真阴性（NT）：正确识别非血管像素为非血管像素的数目。假阳性（PT）：错误识别非血管像素为血管像素的数目。假阴性（NF）：错误识别血管像素为非血管像素的数目。视网膜血管分割结果的定量分析如表2到表6所示，“---”表示相应的比较文献中没有这项指标。表2   HRF_healthy数据库视网膜血管分割结果 HRF_healthy  Method  Sensitivity  FPF  Specificity  Accuracy  15 个数据 文献[8]  0.813947  0.045180  0.954753  0.941633 文献[14]  0.783967  0.028333  0.971667  0.952933 文献[30]  0.763820  0.026307  0.973693  0.954267 文献[36]  0.758573  0.022600  0.977340  0.956787 本文方法  0.794547  0.021100  0.978900  0.961820 表3   HRF_diabetic数据库视网膜血管分割结果 HRF_diabetic  Method  Sensitivity  FPF  Specificity  Accuracy  15 个数据 文献[8]  0.786927  0.067780  0.932227  0.922333 文献[14]  0.742700  0.045133  0.954867  0.941100 文献[30]  0.758940  0.048432  0.951380  0.938160 文献[36]  0.800623  0.057554  0.942446  0.932154 12  计  算  机  学  报                                                                                              2016年 本文方法  0.810653  0.048207  0.951793  0.942147 表4   HRF_glaucoma数据库视网膜血管分割结果 HRF_glaucoma  Method  Sensitivity  FPF  Specificity  Accuracy  15  个数据 文献[8]  0.796627  0.051980  0.948020  0.937867 文献[14]  0.818133  0.051033  0.949033  0.940100 文献[30]  0.753787  0.034707  0.965293  0.950960 文献[36]  0.794800  0.040667  0.959513  0.947940 本文方法  0.811773  0.033053  0.966947  0.956413 表2到表6给出了不同算法在HRF和STARE、DRIVE数据库上不同视网膜的血管分割准确率。在HRF数据库上，从表2、3、4可知本文算法在保证较高敏感度和较低的假阳性率前提下，准确率高于文献[8] 、文献[14]、文献[30]、文献[36]。对于表2HRF_healthy数据，本文算法血管分割的敏感度、准确率分别达到79.4547%、96.1820%；对于表3HRF_diabetic数据，本文算法血管分割的敏感度、准确率分别达到81.0653%、94.2147%；对于表4HRF_glaucoma数据，本文算法血管分割的敏感度、准确率分别达到81.1733%、95.6413%。在STARE数据库上，从表5可知文献[14]的准确率虽然比本文算法、文献[8,9]、文献[15]、  文献[19,20]  、  文献[22]、文献[30]和文献[36]略高，但本文算法有最高的敏感度，且比文献[14]高了8.5661个百分点，敏感度表示成功分割血管像素为血管像素的比率，因此文献[14]存在对血管分割不足。对于STARE数据本文算法的血管分割的敏感度、准确率分别达到79.0860%、95.0340%。在DRIVE数据库上，从表6 可知本文算法明显优于文献[8,9]、文献[14,15]、文献[21]、文献[30]和文献[36]，具有最高的准确率，且敏感度和特异性与文献[20]差距不大，假阳性率也与文献[14]相近，血管分割的敏感度、准确率分别达到75.3535%、95.3565%。因此从表2到表6血管分割指标的定量展示可知，在HRF和STARE、DRIVE数据库上本文算法优于其它血管分割算法，另外本文算法可实现血管全自动分割，可见具有较好的临床使用价值。表5  STARE数据库视网膜血管分割结果 STAREDatabase  Method  Sensitivity  FPF  Specificity  Accuracy   20 个数据 文献[8]  0.782874  0.056711  0.943289  0.931863 文献[9]  0.707100  ---  0.972900  0.938800 文献[14]  0.705145  0.029033  0.970965  0.952390 文献[15]  ---  ---  ---  0.945200 文献[19]  ---  ---  ---  0.952100 文献[20]  0.754400  ---  0.981950  0.952100 文献[22]  0.718700  ---  0.976700  0.950900 文献[30]  0.758881  0.045363  0.954625  0.939256 文献[36]  0.784321  0.036274  0.963726  0.949800 本文方法  0.790860  0.037025  0.962980  0.950340 表6  DRIVE数据库视网膜血管分割结果 DRIVEDatabase  Method    Sensitivity  FPF  Specificity  Accuracy   20 个数据 文献[8]  0.753585  0.034605  0.965175  0.947020 文献[9]  0.731500  ---  0.972300  0.947600 文献[14]  0.723990  0.025780  0.974215  0.953275 文献[15]  ---  ---  ---  0.952600 文献[19]  ---  ---  ---  0.946100 文献[20]  0.784150  ---  0.979900  0.951200 文献[21]  0.732700  ---  0.971500  0.937550 文献[22]  0.735400  ---  0.978900  0.947700 文献[30]  0.710675  0.031215  0.968785  0.946630 文献[36]  0.734815  0.026720  0.973280  0.952765 本文方法  0.753535  0.027620  0.972380  0.953565 图11给出了本文算法和其它算法的F-measure性能评价曲线，可从(a)、(b)、(c)看出，在HRF数据库上，不论视网膜是否发生病变(糖尿病和青光眼)，本文的F-measure曲线优越于文献[8]、文献[14] 、文献[30]、文献[36]。对于HRF_healthy数据，其F-measure平均值分别为0.7954、0.7283、0.7508、0.7591、0.7659；对于HRF_diabetic数据，其F-measure平均值分别为0.6664、0.6052、0.6379、0.6456、0.6006；对于HRF_glaucoma数据，其F-measure平均值分别为0.7234、0.6489、0.66613、0.6861、0.6826。从(d)、(e)可以看出在STARE、论文在线出版号  No.173                          梁礼明等：融合形状先验的水平集眼底图像血管分割                                    11  DRIVE数据库上，本文的F-measure曲线也优于文献[8]、文献[14] 、文献[30]、文献[36]。在STARE数据库上，其F-measure平均值分别为0.7251、0.6735、0.7114、0.6679、0.7048。在DRIVE数据库上，其F-measure平均值分别为0.7377、0.7136、0.7339、0.7015、0.7271。另外，从图(b)和(d)知，五种血管分割算法的F-measure波动都较大，这是因为15个HRF_diabetic数据全是由糖尿病的视网膜数据，20个STARE数据中有10个病变视网膜数据，两种数据库的病灶较多，致使分割算法的性能都受到一定影响。其中STARE数据的im0001、im0002、  im0044视网膜具有较严重的硬性渗出、im0003具有较严重的软性渗出。1 2 3 4 5 6 7 8 9 10 11 12 13 14 150.20.250.30.350.40.450.50.550.60.650.70.750.80.850.90.951视网膜数据F-measure  本文算法文献[8]文献[36]文献[30]文献[14]       1 2 3 4 5 6 7 8 9 10 11 12 13 14 150.20.250.30.350.40.450.50.550.60.650.70.750.80.850.90.951视网膜数据F-measure  本文算法文献[8]文献[36]文献[30]文献[14]                             (a)HRF_healthy数据                                                                (b) HRF_diabetic数据 1 2 3 4 5 6 7 8 9 10 11 12 13 14 150.20.250.30.350.40.450.50.550.60.650.70.750.80.850.90.951视网膜数据F-measure  本文算法文献[8]文献[36]文献[30]文献[14]       1 2 3 4 5 6 7 8 9 10 11 12 13 14 150.20.250.30.350.40.450.50.550.60.650.70.750.80.850.90.951视网膜数据F-measure  本文算法文献[8]文献[36]文献[30]文献[14] (c)HRF_glaucoma数据                                                                     (d)STARE数据 1 2 3 4 5 6 7 8 9 10 11 12 13 14 150.20.250.30.350.40.450.50.550.60.650.70.750.80.850.90.951视网膜数据F-measure  本文算法文献[8]文献[36]文献[30]文献[14]                                                                                    (e) DRIVE数据 图11   不同算法的F-measure性能评价曲线 6 结论 视网膜血管的定量和定性分析在医学临床上具重大的实际价值。但经典的Hessian矩阵[6]存在血管分割尺寸过宽、血管交叉处和相邻血管难分割等问题，对此本文提出了融合形状先验的水平集眼底图像血管分割算法。首先利用Hessian矩阵对血管进行粗分割；然后截取不同阈值下的粗分割结果初始化水平集函数、构建能量模型的形状约束项和面积约束项的边缘函数；最后成功构建出包含区域能量拟合信息和形状先验的视网膜血管分14  计  算  机  学  报                                                                                              2016年 割模型。本文也利用视网膜粗略的背景方差设置形状约束项、面积约束项、长度惩罚项的权重系数，进一步克服图像噪声和固定时间开销，并使用连通域面积和宽、高信息构建相应算子去除伪影与病灶。本文视网膜血管分割算法分别在HRF和STARE、DRIVE数据库上进行了验证，其中在HRF_healthy、HRF_diabetic、HRF_glaucoma数据上准确率分别达到96.182%、94.2147%、95.6413%；在STARE、DRIVE数据库上(20个视网膜数据，包含正常人、糖尿病、青光眼)准确率达到95.034%、95.3565%。从实验结果的定性分析和定量分析可知，本文的视网膜血管分割算法优越其它几种算法，因此本文算法在血管分割上具有较强的临床应用价值。 参 考 文 献 [1] Aguiree F, Brown A, Cho N H, et al. IDF Diabetes Atlas. sixth edition. International Diabetes Federation, Basel, Switzerland ,2013. [2] Kirbas C, Quek F. A review of vessel extraction techniques and algorithms. Association for Computing Machinery Computing Surveys, 2004, 36(2): 81-121. [3] Wang  Shuang,  Xu  Liang,  Li  Jian-Jun. Epidemiological  study  on  the  relationship between  retinal  microvascular  abnormalities  and  cardiovascular  disease.  Section Ophthalmol Foreign Med ical Science, 2005, 29(3):145-148. (王爽,  徐亮,  李建军.  视网膜微血管异常与心脑血管疾病关系的流行病学研究.  国外医学:眼科学分册, 2005, 29(3):145-148) [4] Vermeer K A, Van der Schoot J, Lemij H G, et al. Automated segmentation by pixel classification of retinal layers in ophthalmic OCT images. Biomedical optics express, 2011, 2(6): 1743-1756. [5] Winder R J, Morrow P J, McRitchie I N, et al. Algorithms for digital image processing in diabetic retinopathy. Computerized Medical Imaging and Graphics, 2009, 33(8): 608-622. [6] Lorenz C, Carlsen I C, Buzug T M, et al. Multi-scale line segmentation with automatic estimation of width, contrast and tangential direction in 2D and 3D medical images. Lecture Notes in Computer Science, 2006, 1205:233-242.. [7] Sato Y, Nakajima S, Atsumi H, et al. 3D Multi-scale line filter for segmentation and visualization  of  curvilinear  structures  in  medical  images.  Lecture  Notes  in  Computer Science, 1998, 2(2):143-168. [8]  Frangi  A  F,  Niessen  W  J,  Vincken  K  L,  et  al.  Multiscale  vessel  enhancement filtering//Proceedings of the International Conference on Medical Image Computing and Computer-Assisted Intervention. Berlin Heidelberg, 1998: 130-137.   [9] Chen Qian-Qing. Retinal Vessel segmentation Based on Hessian Matrix and Multiscale Analysis[Master  Thesis].  Huazhong  University  of  Science  &  Technology,Wuhan,China, 2012. (陈倩清.  基于黑森矩阵和多尺度分析的视网膜血管分割[硕士学位论文]. 华中科技大学,  武汉，2012) [10] Jiang X, Mojon D. Adaptive local thresholding by verification-based multithreshold probing with application to vessel detection in retinal images. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2003, 25(1): 131-137. [11] Walter T, Klein J C. Segmentation of Color Fundus Images of the Human Retina: Detection  of  the  Optic  Disc  and  the  Vascular  Tree  Using  Morphological  Techniques// Proceedings  of the International  Symposium  on Medical  Data  Analysis. Madrid,  Spain, 2001:282-287. [12]  Zana  F,  Klein  J  C.  Segmentation  of  vessel-like  patterns  using  mathematical morphology and curvature evaluation. IEEE transactions on image processing, 2001, 10(7): 1010-1019. [13] Nguyen  U  T  V,  Bhuiyan  A,  Park  L  A  F,  et  al.  An  effective  retinal  blood  vessel segmentation  method  using  multi-scale  line  detection.  Pattern  recognition,  2013,  46(3): 703-715. [14]  Gao  Xiang-Jun.  Retinal  vessel  segmentation  based  on  multi-scale  line  detection. Science Technology and Engineering, 2013, 13(23):6820-6824. (高向军.  基于多尺度线性检测的视网膜血管分割.  科学技术与工程,  2013, 13(23):6820-6824) [15] Marín D, Aquino A, Gegúndez-Arias M E, et al. A new supervised method for blood vessel  segmentation  in  retinal  images  by  using  gray-level  and  moment  invariants-based features. IEEE transactions on medical imaging, 2011, 30(1): 146-158. [16] Nekovei R, Sun Y. Back-propagation network and its configuration for blood vessel detection in angiograms. IEEE Transactions on Neural Networks, 1995, 6(1): 64-72. [17] Mehrotra A, Tripathi S, Singh K K, et al. Blood Vessel Extraction for retinal images using  morphological  operator  and  KCN  clustering//  Proceedings  of  the  2014  IEEE International Advance Computing Conference. Gurgaon, India, 2014: 1142-1146. [18] Rodriguez R, Alarcón T E, Abad J J. Blood vessel segmentation via neural network in histological images. Journal of Intelligent and Robotic Systems, 2003, 36(4): 451-465. [19] Wang Y, Ji G, Lin P, et al. Retinal vessel segmentation using multiwavelet kernels and multiscale hierarchical decomposition. Pattern Recognition, 2013, 46(8): 2117-2133. [20] Guo Z, Lin P, Ji G, et al. Retinal vessel segmentation using a finite element based binary level set method. Inverse Problems and Imaging, 2014, 8(2): 459-473.   [21] Dizdaro  B,  Ataer-Cansizoglu  E,  Kalpathy-Cramer  J,  et  al.  Level  sets  for  retinal vasculature  segmentation  using  seeds  from  ridges  and  edges  from  phase  maps.  IEEE International Workshop Mach Learn Signal Process, 2012, 44(1):1-6.   [22] Zhao Y Q, Wang X H, Wang X F, et al. Retinal vessels segmentation based on level set and region growing. Pattern Recognition, 2014, 47(7): 2437-2446. [23] Ghalehnovi M, Zahedi E, Fatemizadeh E. Integration of spatial fuzzy clustering with level set for segmentation of 2-D angiogram// Proceedings of the 2014 IEEE Conference on Biomedical Engineering and Sciences. Kuala Lumpur, Malaysia, 2014: 309-314. [24] Kass  M,  Witkin  A,  Terzopoulos  D. Snakes:  Active  contour  models.  International journal of computer vision, 1988, 1(4): 321-331. [25] Xu C, Prince J L. Snakes, shapes, and gradient vector flow. IEEE Transactions on image processing, 1998, 7(3): 359-369. [26]  Collewet  C.  Polar  snakes:  A  fast  and  robust  parametric  active  contour  model// Proceedings of the 2009 16th IEEE International Conference on Image Processing. Cairo, Egypt, 2009: 3013-3016. [27] Piovano J, Rousson M, Papadopoulo T. Efficient segmentation of piecewise smooth images//  Proceedings  of  the  International  Conference  on  Scale  Space  and  Variational Methods in Computer Vision. Ischia, Italy, 2007: 709-720. [28] Kichenassamy S, Kumar A, Olver P, et al. Gradient flows and geometric active contour models//Proceedings of IEEE International Conference on Computer Vision. Cambridge, USA, 1995: 810-815. 论文在线出版号  No.173                          梁礼明等：融合形状先验的水平集眼底图像血管分割                                    15  [29] Kimmel R, Amir A, Bruckstein A M. Finding shortest paths on surfaces using level sets propagation. IEEE Transactions on Pattern Analysis and Machine Intelligence, 1995, 17(6): 635-640. [30] Caselles V, Kimmel R, Sapiro G. Geodesic active contours. International journal of computer vision, 1997, 22(1): 61-79. [31] Li C, Xu C, Gui C, et al. Level set evolution without re-initialization: a new variational formulation// Proceedings of the 2005 IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR'05). San Diego, USA, 2005, 1: 430-436. [32] Li C, Kao C Y, Gore J C, et al. Minimization of region-scalable fitting energy for image segmentation. IEEE transactions on image processing, 2008, 17(10): 1940-1949.   [33] Li C, Kao C Y, Gore J C, et al. Implicit active contours driven by local binary fitting energy//  Proceedings  of  the  2007 IEEE  Conference  on  Computer  Vision  and  Pattern Recognition. Minneapolis, USA, 2007: 1-7.  [34] Li C, Wang L, Kao C Y, et al. Brain MR image segmentation by minimizing scalable neighborhood intensity fitting energy: A multiphase level set approach// Proceedings of the International Society for Magnetic Resonance in Medicine. Ontario, Canada, 2008: 556. [35] Vese L A, Chan T F. A multiphase level set framework for image segmentation using the  Mumford  and  Shah  model.  International  journal  of  computer  vision,  2002,  50(3): 271-293. [36] Chan T F, Vese L A. Active contours without edges. IEEE Transactions on image processing, 2001, 10(2): 266-277. [37] Osher S, Sethian J A. Fronts propagating with curvature-dependent speed: algorithms based  on  Hamilton-Jacobi  formulations.  Journal  of  computational  physics,  1988,  79(1): 12-49.                                                                                                                                                  LIANGLi-Ming,born in 1967, M.S. ,Professor, Ph.D.supervisor. His current research interests include medicalimage analysis and computer vision.   HUANG Chao-Lin, born in 1989,M.S.candidate. His research direction isthe partial differential equation method of medical image processing.  SHI Fei, born in 1980, Ph.D.,lecturer. Her  research interests include medical image processing and analysis.  WU Jian, born in 1991,M.S.candidate. His research interests include computervision and medical image processing. JIANG Hong-jiu, born in 1993, undergraduate student. His main  research interest is medical image processing based on level set.   CHENXin-Jian, born in1979,Ph.D.,professor. Hisresearch interests include attern recognitionmachinelearning and medical image processing and analysis .  Background In ophthalmology, retinal vessel lesions canbe caused by many diseases, such as glaucoma, macular, hypertension,arteriosclerosis and diabetes. Retinal vessel images can be achieved by digitalimage processing. They can be analysed qualitatively and quantitatively to helpdoctors diagnose illnesses and make treatment plans. Thus, we can reduce the probabilitiesof blindness and retinal lesions on patients. 2D true-color retinal images canbe fetched by fundus cameras directly. These images have a great value inmedical application, since microvessel and other lesions on the retina can beclearly presented on them. Retinal vessel segmentation methods hasgreatly developed in recent years. However, they all have followingdisadvantages: 1. Fractured microvessel in the result. 2. Segmentation error ofcrossing vessels. 3. Excessively huge size of the vessel in the result. 4. Lessrobustness to the pathological retina.  Based on the current work, we propose aretinal vessel segmentation model using level set combined with shape priori.Hessian matrix is used to complete the initial segmentation, to initialize thelevel set function and to build the shape priori. The evolution of the levelset function can be driven by the data term to overcome the segmentation errorof the crossing vessels. Finally, geometry operators based on the connecteddomain are used to eliminate artifacts and lesions. As the experiment resultsshown, our method achieves a better segmentation result compared with otherexisting methods. This work was supported in part by theNational Basic Research Program14  计  算  机  学  报                                                                                              2016年  of China (973 Program) under Grant2014CB748600, and in part by the National Natural Science Foundation of China(NSFC) under Grant 81371629, Grant 6140293, Grant 61401294, Grant 81401451, andGrant 81401472.