基于眼底视网膜血管网络模型的血管提取方法

第 33 卷第 4 期 1999 年 7 月 浙　　江　　大　　学　　学　　报 Jou rnal of Zhejiang U n iversity (工学版) (Engineering Science) V o l. 33 №. 4 Ju l. 1999 收稿日期: 1996211219 基金项目: 浙江省自然科学基金资助项目 (JD 9503) 作者简介: 许　雷 (1970- ) , 男, 浙江杭州人, 浙江大学博士, 从事数字图像处理研究Λ 文章编号: 025329861 (1999) 0420386207 基于眼底视网膜血管网络模型 的血管提取方法 许　雷1, 郑筱祥1, 俞　锋2, 姜节凯2 (1. 浙江大学 生命科学与医学工程学系, 浙江 杭州 310027; 2. 浙江大学 第二附属医院, 浙江 杭州 310005) 摘　要: 在分析正常与病理状况下眼底图像特点的基础上, 建立了眼底血管网络模型Λ根据该模型, 首先 使用 SED (signed edge detection)与RO SE (open ing w ith linear ro ta t ing structu ring elem ent)算法提取出 血管的粗骨架, 然后利用该骨架选出使用LO G (Ρ= 1. 4) 算子计算过零点, 并加以修正得到的精确血管 边沿点Λ 结果明: 用本文方法提取出的血管边沿连续性好, 位置定位精度高, 且具有虚假边沿产生少、 抗干扰能力强的优点Λ 关键词: 视网膜血管; 边沿提取; 旋转式扫描 中图分类号: R 318　　　　　文献标识码: A 0　前　言 高血压、脑血管硬化、冠状动脉粥样硬化等心脑血管疾病是目前我国老年人死亡和致残的主要 原因Λ 此类疾病损伤的组织水平首先是在微循环和微血管层次的变化Λ 眼底视网膜微血管是人体 唯一可以非创伤性直接观察的较深层的微血管, 它的改变程度与高血压等疾病的病程, 严重程度及 愈后情况密切相关[1 ]Λ考虑到眼底测量结果的准确性、客观性、可重复性及大批量眼底普查的需要, 对眼底图像的自动定量分析尤其是与高血压密切相关的视网膜血管改变的定量分析提出了迫切要 求Λ 为此需采用有效的方法提取血管进行分析Λ 从直方图选取全局阈值, 然后加以二值化的血管提取方法由于受局部光照不匀, 散焦及噪声等 的影响, 使得分割结果很差Λ 采用局部自适应二值化的方法则受眼底图像在生理上的复杂性的影 响, 也达不到良好的效果Λ 采用梯度法 (Sobel 算子等) 提取边沿时, 对具有理想边沿特性的图像结 果是好的, 但同时对噪声也加以放大, 且数字图像的非连续性得以加强Λ由于血管截面近似为椭圆, 加上经照相、数字化等多因素的影响使得实际获得的图像中的边沿被模糊化, 因此梯度法并不适 用Λ 为抑制噪声, 提高边沿提取的连续性,M arr 等提出了基于ý 2 (G á I ( i, j ) ) 过零点的边沿检测 算法 (LO G 算子) Ζ但是该算法在选用小参数 Ρ时, 虽边沿提取的位置精度较高, 但同时也检测出过 多 的细小变化, 并形成许多不该有的闭合环, 易与圆点状出血斑和微血管瘤相混淆ΖF. U lup inar [2 ] 针对零交叉并不总是对应事实上的边沿, 也并不总是准确指示边沿位置的问 err = ∃kd Ρ2 , © 1994-2009 China Academic Journal Electronic Publishing House. All rights reserved. http://www.cnki.net 对边沿进行了修正ΖCanny J. [3 ] 提出了适用于任意形状边沿提取的最优算子 (以 ûý 2 (G á I ( i, j ) ) û 局部极大值点作为沿点) 的思想, 但由于使用的是高斯滤波器, 依然存在着边沿位置不准的问 题Λ利用形态学梯度算子[4 ]则具有方法简单, 速度快的优点, 但只适以只含椒盐噪声的图像Λ以张弛 法抽取具有线状形态的血管时, 由于动静脉交叉和血管分叉, 往往使得一枝血管被截断成几段Λ 采 用H ueckel 算子[5 ] , 以拟合的方法提取血管时, 具有对噪声不敏感, 在高纹理区亦能提取出边沿的 优点, 但其计算量过于庞大Λ考虑被抽取血管本身所具有的特点[6 ] , 即具有线状, 血管左边的梯度方 向和右边界的梯度方向正好相反, T ascin i G. 等[7 ]提出了基于梯度方向性的边沿提取方法 (SED ) Λ 然而高血压Ê 期后的血管边界是十分模糊的, 其过度低的对比度使得跟踪无法进行Λ为增强血管与 背景的对比度, R angayyan R. M. 等提出了在频域中增强具有线状特性血管的方法Λ但由于眼底渗 出物的存在, 结果也不好ΛChaudu ri S. 等[8 ]采用 12 个不同方向的高斯模板对眼底图像进行匹配滤 波以增强血管部分Λ 但是该方法对高血压Ê 期患者 (中轴反光亢进)的眼底图像不能很好地加以增 强Λ 再加上它使用固定尺寸 (13 p ixel×9 p ixel)的高斯模板, 所以对管径变化显著及血管高度扭曲 的图像效果不佳Λ 若采用不同尺寸的高斯模板将使得计算量骤增Λ 从以上分析可知, 单纯采用 SED (提取出的血管连续性较差, 抗干扰能力强)或LO G (血管管径 的微小变化能够被刻划出, 提取出的血管边沿连续性好, 但过多的非血管边沿细节也被检测出) 方 法均不能得到理想的结果, 需将两者的优点结合起来Λ 另外, 无论是 SED 还是LO G 方法所形成的 非血管边沿点都过多, 不利于血管跟踪时血管参数的自动测量Λ 为此, 本文提出了下述基于模型的 改进方法Λ 1　方法与材料 1. 1　眼底血管网络模型 通过分析正常和病理状况下眼底图像的特点 (正常人的动静脉管径比AV CRN≈ 0. 7; 高血压 性视网膜动脉硬化É 级为AV CR °É = 0. 640±0. 107; 动脉硬化Ê 级为AV CR °Ê = 0. 635±0. 106Λ 动 静脉交叉处 (Salu s 征)静脉锥度: V CR °É = 7. 86±7. 86; V CR °Ê = 18. 76±7. 26, 我们初步建立了一个 眼底血管网络模型, 并由此得到一个完全自动化的眼底血管参数测量方法Λ (1) 所有血管均从视神经乳头出发, 以放射状扩展至眼底周边部, 这为旋转式扫描的血管提取 提供依据 (从视神经乳头出发的血管数在不同个体之间为 4～ 16) Λ 管径变化 (近视神经乳头处管 径: CA = 100 Λm ; CV = 200 Λm ) 是一个渐变的过程, 越到周边部越小 (CA = 8～ 15 Λm ; CV = 8～ 20 Λm ) Ζ (2)视神经乳头直径在不同个体 (青春期后)之间基本保持不变Λ 可以此为标准, 来衡量不同放 大倍数下血管的管径Λ视神经乳头形状近似为椭圆, 平均灰度最高, 可采用中轴变换进行快速定位Λ (3)非荧光图像的血管灰度比周围组织低, 荧光图像的血管灰度比周围组织高Λ (4)所有的血管分叉均为二分叉 (便于以二叉树数据结构存储及血管跟踪与遍历) Λ (5) 正常人血管无渗出, 主干血管迂曲度较小Λ 血管径向灰度曲线近似为两个高斯曲线之差: IV (x ) = I b0 + A exp [ - (x - Λ2) 2ö2Ρ22 ] - K exp [ - (x - Λ1) 2ö2Ρ21 ], I b0 为背景灰度Λ 1. 2　血管提取步骤 血管提取流程简图如图 1 所示, 具体步骤如下: (1) 用激光扫描仪将眼底照相机摄制的眼底照片数字化, 并以W indow s D IB (640 p ixel×480 p ixel, 256 级灰度)位图文件存入计算机Λ (2)对数字化后的原始图像进行平滑滤波以消除高斯加性或乘性噪声Λ 783第 4 期　　　　　　　　许　雷等: 基于眼底视网膜血管网络模型的血管提取方法 © 1994-2009 China Academic Journal Electronic Publishing House. All rights reserved. http://www.cnki.net (3)背景信号去除及图像增强Λ 因血管具有良好的形态学特征, 适合于基于数学形态学的血管 图像增强[6 ]Λ首先将滤波过的图像 f 对圆形结构元g (r = 20 p ixel) 进行开运算以产生背景图像, 即 O (f , g ) = (f ( g ) Ý g , 其中 (f Ý g ) (x ) = m axu∈Dx - ug∈D f {f (x - u ) + g (u ) }; (f ( g ) (x ) = m inu∈Dx + ug∈D f {f (x + u ) - g (u ) }; f : D f < Z 2; R f < Z ; g. D g < Z 2Ζ然后对减去背景后的图像进行基于 视觉非线性特性的自适应图像增强 (JGA CE) [9 ] , 得增强后的图像为: I 0 = JGA CE (f - O (f , g ) ) Ζ 利用该图像增强方法可避免全局增强方法引起的局部对比度过提升或欠提升以及平滑区噪声过度 放大的问题Ζ 图 1　 血管提取流程简图 (4) 以视神经乳头处为中心, 建立标系, 如图 2 所示Ζ采用圆弧的DDA 算法[10 ] 生成圆弧扫描 线Λ对第一象限圆弧找描线上的每一点, 采用 Sobel0, 6, 7 号模块 (图 2) 分别对其求卷积, 选取最大 输出幅度值作为该点的边沿值Λ因血管灰度小于背景灰度, 所以, 顺时针扫描时血管左边界产生负 峰值, 右边界产生正峰值 (图 3a) Λ同理, 对第二象限圆弧扫描线上的每一点分别使用 0, 1, 2 号模板 对其求卷积Λ第三象限使用 2, 3, 4 号模板, 第四象限使用 4, 5, 6 号模板Λ经先负后正的模板 (宽 13 p ixel) 匹配后正峰值表示血管, 负峰值表示血管之间的背景 (图 3b) Ζ该方法能显著降低血管错误划 分的概率 (图 2 所示: r = 75 p ixel 圆上, 明亮部分表示血管) Ζ (5) 对经匹配的血管图像进行二值化Ζ结果除血管部分外, 网膜部分也产生许多短小的非血管 片段Ζ采用二值化的血管图像与 8 个方向 (0°, 26. 6°, 45°, 63. 4°, 90°, 116. 6°, 135°, 153. 4°) 线形结构 元 (长为 18 p ixel) 所组成的旋转结构元RO SE 进行开运算[11 ] , 目的在于保留长而直的主干血管部 分, 去除短且弯曲的非血管部分Λ由此得到去除干扰的主干血管粗骨架Λ在动静脉交叉处, 血管往往 间断 (这既可能是由于血管对比度差引起, 也可能是由于高血压 Ê 期患者静脉隐匿引起 ——Salu s 征) Λ因此我们采用膨胀算子对二值化的图像进行膨胀, 由此得到一个连通性较好的血管粗骨架 (图 4e) Λ然后对该血管粗骨架图进行十字形中值滤波, 以平滑血管边界Λ接着采用安全点细化算法 SPTA [12 ] 血管中心线的提取 (图 4f) Λ (6) 采用血管的粗骨架与采用M arr2H ildreth 算子 (LO G) 计算过零点得到的边沿 (边沿点及边 沿位置修正, 参阅文献[ 2 ]) 相“与”的方法, 从而得到既具有较好的血管连通性, 边沿位置定位精度 高, 又不至于有过多细节形成的血管边沿图Λ同时将该边沿图与 Sgn (SED ( I ( i, j ) ) ) 相乘以保留血 管边沿的符号信息, 以便于血管跟踪时血管参数计算及降低错误的概率Ζ 883 浙　江　大　学　学　报 (工学版)　　　　　　　　　1999 年 © 1994-2009 China Academic Journal Electronic Publishing House. All rights reserved. http://www.cnki.net 图 2　 旋转式扫描血管提取示意图 a. r = 75 p ixel 圆上各点与 Sobel 模板卷积后的输出; 　　b. 经过模板匹配后的输出 图 3　r = 75 p ixel 圆上 SED ( I ( i, j ) ) á Tem 的灰度曲线图 (7) 血管跟踪采用跟踪主干血管中心线, 同时在当前血管中心线走向的垂直方向及垂直方向 的附近 (± 45°, 血管粗骨架内) 寻找连接左边沿点、当前跟踪点、右边沿点这三点距离最短的连线 (若不能找到基于LO G 算子过零点的边沿点, 则以 SED 算法得到的局部梯度幅度最大值点作为边 沿点[13 ]) Λ以该连线长度作为当前跟踪点处血管的管径, 对测得的二级分支后的动脉及静脉主干管 径求平均值, 两者之比就是动静脉管径比Λ类似可得其它血管参数 (Gunn 征, 血管迂曲度等) 自动 测量的结果Λ上述血管提取用公式表示Λ ① 血管边沿点为: LR ( i, j ) = Ro se (D ila t ion (B in (SED ( I ( i, j ) ) á T em ) ) ) õ Sgn (SED ( I ( i, j ) ) ) õ (LO G ( I ( i, j ) ) ) ; ②血管骨架图为: T h in ( i, j ) = SPTA (M F ilter (D ila t ion (B in (SED ( I ( i, j ) ) á T em ) ) ) ). 其中 I ( i, j ) 为输入图像; LR ( i, j ) > 0代表右边沿点,L R ( i, j ) < 0代表左边沿点 (顺时针扫描, 非荧光图 像) ;M F ilter () 表示十字型中值滤波; D ila t ion () 为膨胀算子; B in () 为二值化算子; Sgn () 为符号算 子; T em 表示 SED 算法中先负后正的模板; á 表示卷积Ζ 983第 4 期　　　　　　　　许　雷等: 基于眼底视网膜血管网络模型的血管提取方法 © 1994-2009 China Academic Journal Electronic Publishing House. All rights reserved. http://www.cnki.net 2　结果与讨论 为了验证本文方法的有效性, 我们在 PEN T IUM 2S586 型计算机 (120M 主频, 16M 内存) 上进 行了大量的对比实验Ζ结果如下, 血管的错误划分率 (错误划分血管像素数 ö血管部分像素数) : 全 局二值化方法为E G = 427. 4◊ ; 局部二值化法E L = 516. 8◊ ; 梯度法E Sobel = 160. 3◊ ; LO G 法E LO G = 209. 6◊ ; 本文方法 E SLR = 5. 577◊ Ζ边沿宽度:W Sobel = 1～ 3 p ixel; 形态学方法W mo rph = 2～ 5 p ixel;W SLR = 1 p ixelΖ处理时间: 文献[ 8 ] 法 T Gaussian = 10. 7 m in; T SLR = 4. 3 m in; T LO G = 2. 3 m in; T Sobel = 0. 73 m inΖ图 4 为对照图像Ζ其中: 图 4a 是原始图像 (440 p ixel × 420 p ixel × 8 b its) , 其背 景灰度高度不均, 血管对比度很差Ζ图 4b 是采用O tsu 法确定全局阈值 (T = 89) , 然后加以二值化 的结果Ζ由图4b 可见, 大量的非血管部分因背景不均而被错误地分割成血管部分Ζ图4c是采用二尺 度 (W 1 = 20,W 2 = 6) 局部自适应二值化的血管提取结果Ζ图 4d 是采用 Sobel 算子边沿提取的结 果Λ原始图像的噪声被严重放大Λ图 4e 是经 SED、模板匹配、二值化, 膨胀一次再经RO SE 算法 (结 构元长度为 18 p ixel) 处理后的结果Λ显然, 非血管部分得到了显著去除Λ图 4g 是采用LO G 算子 (17 p ixel × 17 p ixel, Ρ = 1. 4) 提取边沿, 并加以修正[2 ] 的结果Ζ显然在如此多的细节被提取出的图像 中是难以进行血管跟踪的Ζ图4h 是本文的提出方法的结果Λ图4i是图4h 经血管跟踪并与图4a 叠加 后的结果Ζ由图可见: 血管边界为单像素宽, 无需细化处理Ζ主干血管边沿连续性好, 非常有利于血 管跟踪的进行Ζ非血管边沿得到了大量清除(其中, 传统方法难以消除的斜纹干扰得到彻底清除, 这主要 是因为该方法充分运用了眼底血管网络模型本身所具有的特点, 采用旋转扫描血管提取的结果) Ζ a. 原始图像 b. 全局二值化 c. 局部自适应二值化 d. 梯度法边沿提取 e. 血管粗骨架 f. 血管中心线 093 浙　江　大　学　学　报 (工学版)　　　　　　　　　1999 年 © 1994-2009 China Academic Journal Electronic Publishing House. All rights reserved. http://www.cnki.net g. LO G 算子边沿提取 h. 本文方法的边沿提取 i. 结果在原图上叠加 图 4　对照图像 从眼底血管网络本身所具有的特点出发, 采用以SED 算法与RO SE 算法提取出的血管粗骨架 去选择使用LO G 算子得到的精确血管边沿点, 从而得到抗干扰能力强, 边沿连续性好, 血管边沿 定位精度高的处理结果Λ对照实验的结果亦表明: 本文方法能够处理传统方法不适用的低质量眼底 图像, 并具有不受血管中央反光和附近平行血管及血管交叉影响的优点Λ该方法为视网膜血管网络 参数的自动测量及图像拼接提供了良好的先决条件Λ 参考文献: [ 1 ]　张惠蓉Λ眼微循环及其相关疾病 [M ]Λ北京: 北京医科大学, 中国协和医科大学联合出版社, 1993Λ164～ 170. [2 ]　U lup inar F,M edion i G. R efin ing edges detected by L oG operato r [J ]. Comput V is Graph and Im age P rocess, 1990, 51: 275～ 298. [3 ]　Canny J. A computational app roach to edge detection [J ]. IEEE T rans, PAM I, 1986, 8: 679～ 698. [4 ] 　Peng J , R usch Ph. M o rpho logical filters and edge detection app licat ion to m edical im aging [J ]. A nnual In ternational Conference of the IEEE Engineering In M edicine and B io logy Society, 1991, 13 (1) : 0251～ 0252. [5 ]　荆仁杰, 叶秀清, 徐胜荣, 等Λ 计算机图像处理 [M ]. 杭州: 浙江大学出版社, 1990Λ 245～ 261. [6 ]　H uang C C, L i C C, Fan N , et a l. A fast mo rpho logical filter fo r enhancem ent of angiograph ic im ages [J ]. A nnual In ternational Conference of the IEEE Engineering in M edicine and B io logy Society, 1991, 13 (1) : 0229 ～ 0230. [7 ]　T ascin i G, Passerin i G, Pu lit i P, et a l . R etina vascu lar netw o rk recogn it ion [J ]. P roc SP IE, 1993, 1898: 322～ 329. [ 8 ]　Chauduri S, Chatterjee S, Katz N , et a l . D etection of b lood vessels in ret inal im ages using tw o2dim ensional m atched filters [J ]. IEEE T rans M ed Im aging, 1989, 8: 263～ 269. [ 9 ]　T 2L J i, SundareshanM K, Roeh rig H. A dap tive im age con trast enhancem ent based on hum an visual p ropert ies [J ]. IEEE T ransM ed Im aging, 1994, 13: 573～ 586. [10 ]　金廷赞. 计算机图形学 [M ]. 杭州: 浙江大学出版社, 1988. 137～ 138. [ 11 ]　B rett D T hack ray,A lan C N elson. Sem i2au tom atic segm entat ion of vascu lar netw o rk im ages using a ro ta t ing structu ring elem ent (RO SE) w ith m athem atical mo rpho logy and dual featu re th resho lding [J ]. IEEE T rans M ed Im aging, 1993, 12: 385～ 392. [12 ]　N accache N J , Sh inghal R. SPTA : A p ropo sed algo rithm fo r th inn ing b inary patterns [J ]. IEEE T rans System ,M an, and cybernetics, 1984, SM C214 (3) : 409～ 418. [ 13 ]　A ggarw al S J , Y ip C Y J , D iller K R , et a l . T rack ing of vessel diam eter fluctuations using digita l im age analysis [J ]. B iom edical Im age P rocessing, SP IE, 1990, 1245: 182～ 191. 193第 4 期　　　　　　　　许　雷等: 基于眼底视网膜血管网络模型的血管提取方法 © 1994-2009 China Academic Journal Electronic Publishing House. All rights reserved. http://www.cnki.net Research on retina l vessels detection m ethod based on the m odel of the retina vascular network XU L ei1, ZH EN G X iao2x iang1, YU Feng2, J IAN G J ie2kai2 (1. D ep t. of L ife Science and B iom edical Engineering, Zhejiang U niv. , H angzhou 310027, Ch ina; 2. T he Second B ranch Ho sp ita l of Zhejiang U niv. , H angzhou 310027, Ch ina) Abstract: T he p rogression ra te and severity of patho logica l changes of the ret ina l vessels can be u sed in the m ass screen ing of card iocereb ra l vascu lar d iseases. In th is paper, w e p resen t a novel a lgo rithm that base on the model of the ret ina vascu lar netw o rk and resu lt in excellen t detect ion of ret ina l vessels by com b inat ion of SED (signed edge detect ion, fo r labelling lef t and righ t edges and fo r locat ing the“best”zero2cro ssings (ZC s) of vessel boundary) a lgo rithm , RO SE (open ing w ith linear ro ta t ing structu ring elem en t, fo r elim inat ing nonvessel segm en ts) a lgo rithm , and LO G ( lap lacian of gau ssian opera to r, fo r removing no ise and p roducing the ZC s fo rm the im age ) opera to r. A deta iled perfo rm ance evaluat ion together w ith the ex ist ing techn iques is g iven to demon stra te the st rong featu res of ou r m ethod. Key words: ret ina l vessels; edge detect ion; ro tary scann ing (责任编辑: 陈　波) 下期发表摘要预告 交互式有线电视系统机顶盒的研制 陈国玺, 周　峰, 唐慧明, 姚庆栋 (浙江大学 信息与电子工程学系, 浙江 杭州　310027) 摘　要: 提出了交互式有线电视系统机顶盒的设计, 机顶盒具备交互道和数字电视通道, 符合 M PEG2 数字电视标准, 并针对多媒体通讯的特点Λ采用M HN P 模型, 开发了超媒体文档和用户端的浏览 器, 机顶盒是交互式有线电视系统中的一项关键技术, 它的研制必然给家用电视带来一场革命Λ 关键词: 交互式电视; 机顶盒; 多媒体; 超媒体 293 浙　江　大　学　学　报 (工学版)　　　　　　　　　1999 年 © 1994-2009 China Academic Journal Electronic Publishing House. All rights reserved. http://www.cnki.net