一种融合形状先验的水平集视网膜血管图像分割方法
【专利摘要】本发明涉及一种融合形状先验的水平集视网膜血管图像分割方法，包括：(1)利用形态学算子和高斯卷积增强视网膜血管图像；(2)采用Hessian矩阵的各向异性特性和改进的血管响应函数粗略分割视网膜血管图像，并作为形状约束和初始化信息；(3)使用形状先验和视网膜图像数据信息构建一个包含局部区域能量拟合项、形状约束项、水平集函数正则性维持项、长度惩罚项、加权面积约束项的视网膜血管分割的水平集模型。本发明的分割结果具有较高的准确性，能够替代手动分割，对于临床相关眼科疾病的诊断与治疗能起到重要的辅助作用和具有较强的临床应用价值。
【专利说明】
一种融合形状先验的水平集视网膜血管图像分割方法
技术领域
[0001] 本发明涉及一种基于水平集模型的视网膜血管图像分割方法，解决了现有模型存 在相邻血管易相连、血管过宽、细小血管易断裂、血管交叉处分割不足等问题。
【背景技术】
[0002] 视网膜是脑部神经组织的延伸，具有复杂的多层次组织结构，其血管病变是致盲 的重要原因之一。水平集方法是解决曲线演化问题的一种强有力的工具，其拓扑适应性强。 它能提供快速的、高准确率的视网膜血管提取方法，为临床眼科医生对疾病的诊断和治疗 提供帮助。在眼科学领域，视网膜血管数量、分支、角度、宽度等信息均可作为与视网膜血管 相关疾病的诊断依据，这为利用数字图像处理技术分割视网膜血管，定性和定量分析判断 患者病情以及研究病理提供了基础。然而目前眼科医生基本上是采用手动方式对视网膜病 变进行定量分析，主观性强，无法保证准确性和一致性。
[0003] 目前血管分割方法多种，如多尺度血管增强滤波、多阈值的血管检测、以形态学为 基础技术血管分割、使用神经网络的血管分割算法、多尺度层分解和局部自适应阈值血管 分割方法以及基于活动轮廓模型的血管分割等。这些方法存在以下的缺陷：（1)易受高斯卷 积算子的影响，分割结果存在相邻血管易相连、血管过宽、细小血管易断裂、血管交叉处分 割不足等问题，导致用于复杂血管结构的视网膜血管分割时精度较低；（2)存在对图像噪声 过于敏感和难以解决目标与背景灰度值交叉等问题。
【发明内容】
[0004] 本发明的目的是针对现有视网膜血管分割方法的不足，提供一种融合形状先验的 水平集视网膜血管图像分割方法，本发明能有效地解决相邻血管易相连、血管过宽、细小血 管易断裂、血管交叉处分割不足以及对图像噪声过于敏感、目标与背景灰度值交叉等问题。
[0005] -种融合形状先验的水平集视网膜血管图像分割方法，包括三个步骤：
[0006] 步骤1，视网膜血管图像预处理:在对视网膜做边缘扩展基础上采用高斯卷积估计 视网膜背景，以消除视网膜原始边界处的灰度突变，获得了一个背景更加均匀的血管增强 图像；
[0007] 步骤2,视网膜血管图像粗略分割:利用改进的血管响应函数和阈值得到视网膜血 管的粗分割图像，并以此信息初始化水平集函数和构建形状约束项，以克服水平集模型对 初始化敏感以及区域能量拟合项对噪声敏感的问题；
[0008] 步骤3,视网膜血管图像精细分割:应用连通域面积Area和宽W、高Η构建几何算子， 消除连通域面积较小的伪影和病灶，获得高准确率的视网膜血管分割图像。
[0009] 所述的步骤1包括以下三个子步骤：
[0010] (a)选取视网膜图像的绿色通道图像I，利用测地线活动轮廓(Geodesic active contours，GAC)模型自动获得视网膜的"掩模"；
[0011] (b)在获取的"掩模"的基础上对视网膜图像作边缘扩展，获得图像iExten;就是采用 形态学算子，对图像作开运算，消除血管中心亮线，从而获得消除了血管中心亮线的视网膜 边缘扩展图像iExten;
[0012] (C)构建高斯模板，所述的高斯模板尺寸为原始图像尺寸除以14后取整数，小数点 后四舍五入;对图像I Exten做卷积运算，估计视网膜图像背景IBackg，并与IExten做减法运算和 与"掩模"做点乘运算，从而获得增强的视网膜血管图像然后把图像的灰度值 线性拉伸到[0-255]，获得图像1 &_。
[0013] 所述的步骤2包括以下两个子步骤：
[0014] (a)构建血管响应函数，根据Hessian矩阵对图像不同点（x，y)的特征值描述，构建 用以区别血管与背景的描述算子;构建的视网膜血管描述算子J和血管响应函数vd^y)分 别定义为：
[0015] J=(|A1|-|A2|)2(|A1| + |A2|)2 (1)
[0017] 式（1)、式⑵中为二维Hessian矩阵的特征值;β为调整血 管明暗的尺寸参数，σ是高斯核函数尺寸，Jmax是J的最大值；
[0018]由于视网膜血管图像尺寸不一，需在多尺度下计算血管响应函数V(3(X，y)，最后取 各尺度下的最大响应值，其定义如下：
[0019 ] ^(.v, .v) = maxv^.(.v. V) (3)
[0020] 式(3)中〇mir^P〇max为感兴趣的眼底视网膜血管的最小和最大尺寸；
[0021] (b)二值化处理，利用式(2)和式(3)计算v(x，y)，获得血管响应图像IRest^进行二 值化处理;先选取较小的阈值^获得二值化图像，并作为水平集函数φ的初始化;然后 选取略大的阈值(： 2获得二值化图像4IliW，并使用2X2的模板对/^"作开运算和删除连通 域面积小于100的伪影，获得二值化图像Ιβ去构建形状约束项。
[0022] 所述的步骤3包括以下两个子步骤：
[0023] (a)构建视网膜图像血管分割的水平集模型，利用粗略去除血管后视网膜背景方 差设置形状约束项、面积约束项、长度惩罚项的权重系数，进一步克服图像噪声和固定时间 开销;视网膜图像血管分割的水平集模型Ε( Φ，fi，f2)定义如下：
[0024] Ε( Φ ,fi,f2)=ER( Φ ,fi,f2)+nEs( Φ )+uPr( Φ )+vLp( Φ )+yAr( Φ ) (4)
[0025] 式(4)右边分别为局部区域能量拟合项、形状约束项、水平集函数正则性维持项、 长度惩罚项、加权面积约束项，η,μ,ν, γ为权重系数；
[0026]①区域能量拟合项Er( Φ，&彳2)定义为：
[0027] ⑴
[0028] 式(5)中kdx-y)为高斯核函数，〇是高斯核尺寸参数，用以控制局部区域灰度值拟 合范围；f1( X)，i = l，2为局部区域灰度值拟合函数，通过标准梯度下降流法求解能量函数E (Φ，fl，f2)得到：
[0030]式（6)0i，i = 1，2为区域能量拟合项权重系数；Φ为Lipschitz水平集函数， <(办/ = 1，2为区域函数，表示曲线内部和外部，定义ΜΓ以) = Ag⑷=1- Ηε( Φ )为具有正则性的Heaviside函数，其导数为Dirac函数δε( φ );
[0033] 式(7)、式(8)中ε为尺寸参数；
[0034] 由于/^^作为水平集函数的初始化信息，即Μχ,.ν,Ο) = 2卜-24J ;
[0035]②形状约束项Ε5(Φ)能量表达式定义为：
[0036] .⑷=⑷-//..?|: Λ ,；?
[0037] 相应的水平集函数?φ'ν,0) = |二(10)
[0038] 式(9)中识表示视网膜血管先验形状，即在t*时刻 d{x,C) X e Cwlii,jde
[0039] φ(.\\?) - 〇 xe C (II) rd{x,C) xeC?*
[0040] 式（11)中d(x，C)表示视网膜图像中点x到轮廓曲线C的Euclidean距离，xeCoutside 表示点在轮廓线C外面，X e Cinsicb表示点在轮廓线C内部，X e C表示点在轮廓线C上；
[0041 ]③水平集函数正则性维持项Pr( Φ )定义为：
[0042] = ▽詞-1)3办 （12)
[0043] ④长度惩罚项LP( Φ )定义为：
[0044] Lp [φ) - |ω VHt: (^)| dx (13)
[0045] ⑤加权面积约束项Ar( Φ )定义为：
[0046] Ar( Φ )=/Qg(l-He( Φ ))dx (14)
[0047] 式（12)、式（13)和式（14)中Ω为图像全域，g = l/(l + |VG*/Rea|2)为非负单调递减 函数，G为高斯函数，IRea等于Ib点乘IExten;在非血管的背景区域g=l，面积约束项具有一个 较大的惩罚，而在血管附近g-o,面积约束项只获得一个较小的惩罚；
[0048] (b)移除病灶和伪影，利用连通域面积Area和宽W、高Η信息构建相应算子去除伪影 与病灶；
[0049] 当W/H和W X Η同时满足:0.33〈W/H〈3.0和W X H〈4Area时，则可判定整个连通域为病 灶而非血管；当连通域面积Area小于50时，可判定整个连通域为伪影而非血管，从而得以提 高视网膜血管的分割精度。
[0050] 本发明在HRF和STARE、DRIVE数据库上的实验验证，准确率分别达到96.182%、 95.034 %和95.357 %，优于现有分割方法。同时克服了在相邻血管处、血管交叉处和微血管 处其它方法的不足，使分割出的血管结构最为接近金标准和血管真实尺寸值。
【附图说明】
[0051] 图1为本发明实施例中的视网膜图像的绿色通道图像I。
[0052]图2为本发明实施例中自动获得视网膜的"掩模"。
[0053]图3为本发明实施例中消除了血管中心亮线的视网膜边缘扩展图像IExten。
[0054] 图4为本发明实施例中视网膜图像背景估计IBac；kg。
[0055] 图5为本发明实施例中增强的视网膜血管图像Ienhance。
[0056] 图6为本发明实施例中灰度值线性拉伸后获得的图像ICrayS。
[0057] 图7为本发明实施例中获得血管响应图像IResp。
[0058] 图8为本发明实施例中选取阈值C1获得的二值化图像
[0059] 图9为本发明实施例中使用尺寸为2X2的模板对/Law作开运算和删除连通域面 积小于100的伪影，获得的二值化图像Ib。
[0060] 图10为本发明实施例中的12_h视网膜图像伪影和病灶。
[0061 ]图11为本发明实施例中的12_h移除伪影和病灶后分割图像。
【具体实施方式】
[0062]下面结合【具体实施方式】，进一步阐述本发明。
[0063]实验说明：本发明的应用所涉及实施例数据来自于HRF数据库中第12个正常人 (12_h)的视网膜图像。
[0064]本实施例包括三个步骤:视网膜血管图像预处理、血管图像粗略分割和血管图像 精细分割。
[0065] 具体描述如下：
[0066] 1、视网膜血管图像预处理
[0067] (a)选取视网膜图像的绿色通道图像I，利用测地线活动轮廓(Geodesic active contours，GAC)模型自动获得视网膜的"掩模"，如图2所示。
[0068] (b)利用上一步骤(a)获取的视网膜"掩模"信息，对图1做基于镜像对称的边缘扩 展，边缘扩展的尺寸等于下一步骤(c)中的高斯模板尺寸。然后使用板尺寸为5X5形态学算 子，对图像作开运算，消除血管中心亮线，可防止含有区域能量拟合项的水平集模型产生错 误分割，从而获得消除了血管中心亮线的视网膜边缘扩展图像I Ext_如图3所示。
[0069] (c)构建尺寸为100X100的高斯模板，对图像IExten做卷积运算，估计视网膜图像背 景I Backg，如图4所示，并与之做减法运算和与"掩模"做点乘运算，从而获得增强的视网膜血 管图像Ienhana，如图5所示。然后把图像I enhance的灰度值线性拉伸到[0 - 2 5 5 ]，获得图像 I GrayS，如图6所示。
[0070] 2、视网膜血管图像粗略分割
[0071] 视网膜血管图像预粗略分割主要包括以下两个步骤：
[0072] (a)构建血管响应函数
[0073]根据Hessian矩阵对图像不同点（x，y)的特征值描述，构建用以区别血管与背景的 描述算子。本发明构建的视网膜血管描述算子J和血管响应函数wUd)分别定义为：
[0076] 式（15)、式（16)中，用于描述血管，λ^Ρλ；^别为二维Hessian矩阵的特 征值。β为调整血管明暗的尺寸参数，这里取β=0.5;〇是高斯核函数尺寸，Jmax是J的最大值。 上式（16)能最大化血管像素的响应，使血管响应函数vXx，y) - l,而在非血管背景平坦区域 上使血管响应函数νσ(χ,γ)4〇,即[0,1]。
[0077] 由于视网膜血管图像尺寸不一，需在多尺度下计算血管响应函数V(3(X，y)，最后取 各尺度下的最大响应值，其定义如下：
[0078] ^(.v.v) = max^(A-,.v) ?Π)
[0079] 式（17)中〇mir^P〇max为感兴趣的眼底视网膜血管的最小和最大尺寸。这里〇e[l， 11]，步长取2。
[0080] (b)二值化处理
[0081] 利用式（16)和式（17)计算v(x，y)，获得血管响应图像IResp(如图7)后进行二值化 处理。本发明先选取较小的阈值 C1 = 〇. 〇〇〇〇 1获得二值化图像如图8所示，并作为水平 集函数Φ的初始化。然后选取略大的阈值c2 = 0.0001获得二值化图像，并使用尺寸为 2X2的模板对iLrfl.y作开运算和删除连通域面积小于100的伪影，获得二值化图像Ib，如图9 所示。
[0082] 3、视网膜血管图像精细分割
[0083]视网膜血管图像预精细分割主要包括以下两个步骤：
[0084] (a)构建视网膜图像血管分割的水平集模型
[0085] 本发明利用粗略去除血管后视网膜背景方差设置形状约束项、面积约束项、长度 惩罚项的权重系数，进一步克服图像噪声和固定时间开销。视网膜图像血管分割的水平集 模型Ε( Φ，fl，f2)定义如下：
[0086] Ε( Φ ,fi,f2)=ER( Φ ,fi,f2)+nEs( Φ )+uPr( Φ )+vLp( Φ )+yAr( Φ ) (18)
[0087] 式（18)右边分别为局部区域能量拟合项、形状约束项、水平集函数正则性维持项、 长度惩罚项、加权面积约束项。η,μ,ν, γ为权重系数。
[0088]①区域能量拟合项Er( Φ，&彳2)定义为：
[0089]在(為/，/2)=办俯认卜"ν)|/〇')-,/;(.Υ)「Λ<(<ν))4，)?λ. (19)
[0090]式（19)中kdx-y)为高斯核函数，〇为高斯核尺寸参数，用以控制局部区域灰度值 拟合范围；f i(X)，i = 1，2为局部区域灰度值拟合函数，通过标准梯度下降流法求解能量函 数Ε( Φ，fl，f2)得到
[0092] &4 = 1，2为区域能量拟合项权重系数，这里01 = 02 = 1;巾为1^口8〇11；^2水平集函 数，M,W = 1，2为区域函数，表示曲线内部和外部。定义= ，树(# = 1-矿(斗 Ηε( Φ )为具有正则性的Heaviside函数，其导数为Dirac函数δε( φ )。
[0095] 其中ε为尺寸参数，这里取ε = 1。
[0096]由于 <胃作为水平集函数的初始化信息，即#^,0) = 2(1-2/,1^^
[0097]②形状约束项Ε5(Φ)能量表达式定义为：
[0098] Es (φ) = f〇|//'，； [φ] - ??'： dx (23)
[0099] 相应的水平集函数 #(λ·, )，,0) = p ^-° (24) 1 -2. 1 β U
[0100] 式(23)中识表示视网膜血管先验形状，即在t*时刻 φ-，<τ) X e CoutsUk
[0101 ] φ^χ,?*^ = < 0 xeC (25) -d(x,C) x^Cimille
[0102] 其中d(x，C)表示视网膜图像中点x到轮廓曲线C的Euc 1 i dean距离，x e Couts ide表示 点在轮廓线C外面，x e Cinsicb表示点在轮廓线C内部，x e C表示点在轮廓线C上。
[0103]③平集函数正则性维持项PR( Φ )定义为：
[0105] ④长度惩罚项LP( Φ )定义为：
[0106] LP(>5) = Viy,，；(y)|i/v (27)
[0107] ⑤加权面积约束项ΑΚ(Φ)定义为：
[0108] Ακ(Φ)=/Ωδ(1-Ηε(Φ))?χ (28)
[0109] 式（26)、式（27)和式（28)中Ω为图像全域，§ = + 为非负单调递减 函数，G为高斯函数，IRea等于ΙΒ点乘IExten。在非血管的背景区域g=l，面积约束项具有一个 较大的惩罚，而在血管附近g-〇，面积约束项只获得一个较小的惩罚。
[0110] 在视网膜血管分割实验中，时间步长t = 0.1，正则化项系数μ=1、长度惩罚项系数 v=ms、面积约束项系数y=nlogas、形状约束项系数n=ls2(s2表示粗略去除血管后视网膜 背景方差），本实施例中m、n、a、l分别取5、1、7、1。
[0111] (b)移除病灶和伪影
[0112] 本发明利用连通域面积Area和宽、高(W，H)信息构建相应算子去除伪影与病灶。
[0113] 本实施例构建8连通域模板标记水平集分割结果，从而获得每个连通域的面积 Ar ea和宽W、高財言息。
[0114] 当W/H和WXH同时满足：
[0115] 算子①：0.4〈W/H〈2.5
[0116] 算子②：WXH〈3.5Area
[0117] 则整个连通域判定为病灶(非血管）。
[0118] 算子③:Area〈30,则整个连通域判定为伪影(非血管）。
[0119] 利用算子①、算子②和算子③移除伪影和部分病灶，如图10所示，以提高视网膜血 管的分割精度，如图11所示。
【主权项】
1. 一种融合形状先验的水平集视网膜血管图像分割方法，包括Ξ个步骤： 步骤1，视网膜血管图像预处理:在对视网膜做边缘扩展基础上采用高斯卷积估计视网 膜背景，W消除视网膜原始边界处的灰度突变，获得了一个背景更加均匀的血管增强图像； 步骤2,视网膜血管图像粗略分割:利用改进的血管响应函数和阔值得到视网膜血管的 粗分割图像，并W此信息初始化水平集函数和构建形状约束项，W克服水平集模型对初始 化敏感W及区域能量拟合项对噪声敏感的问题； 步骤3，视网膜血管图像精细分割:应用连通域面积Area和宽W、高Η构建几何算子，消除 连通域面积较小的伪影和病灶，获得高准确率的视网膜血管分割图像。2. 根据权利要求1所述的一种融合形状先验的水平集视网膜血管图像分割方法，其特 征是:所述的步骤1包括W下Ξ个子步骤： (a) 选取视网膜图像的绿色通道图像I，利用测地线活动轮廓模型自动获得视网膜的 "掩模'； (b) 在获取的"掩模"的基础上对视网膜图像作边缘扩展，获得图像lExten;就是采用形态 学算子，对图像作开运算，消除血管中屯、亮线，从而获得消除了血管中屯、亮线的视网膜边缘 扩展图像lExten; (C)构建高斯模板，所述的高斯模板尺寸为原始图像尺寸除W14后取整数，小数点后四 舍五入;对图像lExten做卷积运算，估计视网膜图像背景iBackg，并与lExten做减法运算和与"掩 模"做点乘运算，从而获得增强的视网膜血管图像lenhance，然后把图像lenhance的灰度值线性 拉伸到[0-255]，获得图像L·rayS。3. 根据权利要求1所述的一种融合形状先验的水平集视网膜血管图像分割方法，其特 征是:所述的步骤2包括W下两个子步骤： (a) 构建血管响应函数，根据Hessian矩阵对图像不同点(x，y)的特征值描述，构建用W 区别血管与背景的描述算子;构建的视网膜血管描述算子J和血管响应函数vn(x，y)分别定 乂为：式（1)、式(2)中1/Κρ = λ2/λι，λι和λ2分别为二维Hessian矩阵的特征值;β为调整血管明 暗的尺寸参数，σ是高斯核函数尺寸，Jmax是J的最大值； 由于视网膜血管图像尺寸不一，需在多尺度下计算血管响应函数Vn(X，y)，最后取各尺 度下的最大响应值，其定义如下：式(3)中Omin和Omax为感兴趣的眼底视网膜血管的最小和最大尺寸； (b) 二值化处理，利用式(2)和式(3)计算V(X，y)，获得血管响应图像iResp后进行二值化 处理;先选取较小的阔值C1获得二值化图像瑞，并作为水平集函数Φ的初始化;然后选取 略大的阔值C2获得二值化图像巧1。，。，并使用2X2的模板对巧。wy作开运算和删除连通域面 积小于100的伪影，获得二值化图像Ib去构建形状约束项。4.根据权利要求1所述的一种融合形状先验的水平集视网膜血管图像分割方法，其特 征是:所述的步骤3包括W下两个子步骤： (a)构建视网膜图像血管分割的水平集模型，利用粗略去除血管后视网膜背景方差设 置形状约束项、面积约束项、长度惩罚项的权重系数，进一步克服图像噪声和固定时间开 销;视网膜图像血管分割的水平集模型E( Φ，fl，f2)定义如下： E( Φ ,fi,f2)=ER( Φ ,fi,f2)+nEs( Φ )+μΡκ( Φ )+vLp( Φ )+yAr( Φ ) (4) 式(4)右边分别为局部区域能量拟合项、形状约束项、水平集函数正则性维持项、长度 惩罚项、加权面积约束项，η，μ，V，丫为权重系数； ① 区域能量拟合项Er( Φ，fl，f2)定义为：C5) 式(5)中kn(x-y)为高斯核函数，σ是高斯核尺寸参数，用W控制局部区域灰度值拟合范 围；fi(x)，i = l，2为局部区域灰度值拟合函数，通过标准梯度下降流法求解能量函数Ε(Φ, fl,f2)得到：式(6)βι4 = 1，2为区域能量拟合项权重系数；Φ为Lipschitz水平集函数，W;'·'(約，/ = 1，2 为区域函数，表示曲线内部和外部，定义Μ,ι''(辦=(刮.Λ'/;'(抑=1 -/-Γ'(刮;HE ( Φ )为具有 正则性的化aviside函数，其导数为Dirac函数δΕ( Φ );式(7)、式(8)中ε为尺寸参数； 由于這。。。作为水平集函数的初始化信息，即クl·l义.1'，()) = 2?l--2瑞。"y；); ② 形状约束项Es( Φ )能量表达式定义为：式（11)中d(x,C)表示视网膜图像中点X到轮廓曲线C的Euclidean距离，xeCoutside表示 点在轮廓线C外面，X e Cinside表示点在轮廓线C内部，X e C表示点在轮廓线C上； ③水平集函数正则性维持项Pr( Φ )定义为：式（12)、式（13)和式（14)中Ω为图像全域，留= 1/(?十|VG*/J2)为非负单调递减函数， G为高斯函数，iRea等于Ib点乘lExten;在非血管的背景区域g=l，面积约束项具有一个较大的 惩罚，而在血管附近g^O,面积约束项只获得一个较小的惩罚； (b)移除病灶和伪影，利用连通域面积Area和宽W、高Η信息构建相应算子去除伪影与病 灶； 当W/H和W X Η同时满足:0.33<W/H<3.0和W X H<4Area时，则可判定整个连通域为病灶而 非血管；当连通域面积Area小于50时，可判定整个连通域为伪影而非血管，从而得W提高视 网膜血管的分割精度。
【文档编号】G06T7/00GK106097378SQ201610585839
【公开日】2016年11月9日
【申请日】2016年7月24日 公开号201610585839.5, CN 106097378 A, CN 106097378A, CN 201610585839, CN-A-106097378, CN106097378 A, CN106097378A, CN201610585839, CN201610585839.5
【发明人】梁礼明, 黄朝林, 陈新建, 曾璐, 周发助
【申请人】江西理工大学