CN106097378A - 一种融合形状先验的水平集视网膜血管图像分割方法 - Google Patents

Abstract

本发明涉及一种融合形状先验的水平集视网膜血管图像分割方法，包括：(1)利用形态学算子和高斯卷积增强视网膜血管图像；(2)采用Hessian矩阵的各向异性特性和改进的血管响应函数粗略分割视网膜血管图像，并作为形状约束和初始化信息；(3)使用形状先验和视网膜图像数据信息构建一个包含局部区域能量拟合项、形状约束项、水平集函数正则性维持项、长度惩罚项、加权面积约束项的视网膜血管分割的水平集模型。本发明的分割结果具有较高的准确性，能够替代手动分割，对于临床相关眼科疾病的诊断与治疗能起到重要的辅助作用和具有较强的临床应用价值。

Description

一种融合形状先验的水平集视网膜血管图像分割方法

技术领域

本发明涉及一种基于水平集模型的视网膜血管图像分割方法，解决了现有模型存在相邻血管易相连、血管过宽、细小血管易断裂、血管交叉处分割不足等问题。

背景技术

视网膜是脑部神经组织的延伸，具有复杂的多层次组织结构，其血管病变是致盲的重要原因之一。水平集方法是解决曲线演化问题的一种强有力的工具，其拓扑适应性强。它能提供快速的、高准确率的视网膜血管提取方法，为临床眼科医生对疾病的诊断和治疗提供帮助。在眼科学领域，视网膜血管数量、分支、角度、宽度等信息均可作为与视网膜血管相关疾病的诊断依据，这为利用数字图像处理技术分割视网膜血管，定性和定量分析判断患者病情以及研究病理提供了基础。然而目前眼科医生基本上是采用手动方式对视网膜病变进行定量分析，主观性强，无法保证准确性和一致性。

目前血管分割方法多种，如多尺度血管增强滤波、多阈值的血管检测、以形态学为基础技术血管分割、使用神经网络的血管分割算法、多尺度层分解和局部自适应阈值血管分割方法以及基于活动轮廓模型的血管分割等。这些方法存在以下的缺陷：(1)易受高斯卷积算子的影响，分割结果存在相邻血管易相连、血管过宽、细小血管易断裂、血管交叉处分割不足等问题，导致用于复杂血管结构的视网膜血管分割时精度较低；(2)存在对图像噪声过于敏感和难以解决目标与背景灰度值交叉等问题。

发明内容

本发明的目的是针对现有视网膜血管分割方法的不足，提供一种融合形状先验的水平集视网膜血管图像分割方法，本发明能有效地解决相邻血管易相连、血管过宽、细小血管易断裂、血管交叉处分割不足以及对图像噪声过于敏感、目标与背景灰度值交叉等问题。

一种融合形状先验的水平集视网膜血管图像分割方法，包括三个步骤：

步骤1，视网膜血管图像预处理：在对视网膜做边缘扩展基础上采用高斯卷积估计视网膜背景，以消除视网膜原始边界处的灰度突变，获得了一个背景更加均匀的血管增强图像；

步骤2，视网膜血管图像粗略分割：利用改进的血管响应函数和阈值得到视网膜血管的粗分割图像，并以此信息初始化水平集函数和构建形状约束项，以克服水平集模型对初始化敏感以及区域能量拟合项对噪声敏感的问题；

步骤3，视网膜血管图像精细分割：应用连通域面积Area和宽W、高H构建几何算子，消除连通域面积较小的伪影和病灶，获得高准确率的视网膜血管分割图像。

所述的步骤1包括以下三个子步骤：

(a)选取视网膜图像的绿色通道图像I，利用测地线活动轮廓(Geodesic activecontours，GAC)模型自动获得视网膜的“掩模”；

(b)在获取的“掩模”的基础上对视网膜图像作边缘扩展，获得图像I_Exten；就是采用形态学算子，对图像作开运算，消除血管中心亮线，从而获得消除了血管中心亮线的视网膜边缘扩展图像I_Exten；

(c)构建高斯模板，所述的高斯模板尺寸为原始图像尺寸除以14后取整数，小数点后四舍五入；对图像I_Exten做卷积运算，估计视网膜图像背景I_Backg，并与I_Exten做减法运算和与“掩模”做点乘运算，从而获得增强的视网膜血管图像I_enhance，然后把图像I_enhance的灰度值线性拉伸到[0-255]，获得图像I_GrayS。

所述的步骤2包括以下两个子步骤：

(a)构建血管响应函数，根据Hessian矩阵对图像不同点(x,y)的特征值描述，构建用以区别血管与背景的描述算子；构建的视网膜血管描述算子J和血管响应函数ν_σ(x,y)分别定义为：

J＝(|λ₁|-|λ₂|)²(|λ₁|+|λ₂|)² (1)

$v_{\mathrm{\&sigma;}}(x,y)=\left\{\begin{array}{cc}0& {\mathrm{\&lambda;}}_2<0\\ (1-\exp (-\frac{R_{\mathrm{\&beta;}}^2}{2{\mathrm{\&beta;}}^2}\left)\right)(1-\exp (-\frac{J}{\sqrt{J_{\max }}}\left)\right)& {\mathrm{\&lambda;}}_2>=0\end{array}\right.---\left(2\right)$

式(1)、式(2)中1/R_β＝λ₂/λ₁，λ₁和λ₂分别为二维Hessian矩阵的特征值；β为调整血管明暗的尺寸参数，σ是高斯核函数尺寸，J_max是J的最大值；

由于视网膜血管图像尺寸不一，需在多尺度下计算血管响应函数ν_σ(x,y)，最后取各尺度下的最大响应值，其定义如下：

$v(x,y)=\underset{{\mathrm{\&sigma;}}_{\min }\&le;\mathrm{\&sigma;}\&le;{\mathrm{\&sigma;}}_{max}}{{\mathrm{maxv}}_{\mathrm{\&sigma;}}(x,y)}---\left(3\right)$

式(3)中σ_min和σ_max为感兴趣的眼底视网膜血管的最小和最大尺寸；

(b)二值化处理，利用式(2)和式(3)计算ν(x,y)，获得血管响应图像I_Resp后进行二值化处理；先选取较小的阈值c₁获得二值化图像并作为水平集函数φ的初始化；然后选取略大的阈值c₂获得二值化图像并使用2×2的模板对作开运算和删除连通域面积小于100的伪影，获得二值化图像I_B去构建形状约束项。

所述的步骤3包括以下两个子步骤：

(a)构建视网膜图像血管分割的水平集模型，利用粗略去除血管后视网膜背景方差设置形状约束项、面积约束项、长度惩罚项的权重系数，进一步克服图像噪声和固定时间开销；视网膜图像血管分割的水平集模型E(φ,f₁,f₂)定义如下：

E(φ,f₁,f₂)＝E_R(φ,f₁,f₂)+ηE_S(φ)+μP_R(φ)+νL_P(φ)+γA_R(φ) (4)

式(4)右边分别为局部区域能量拟合项、形状约束项、水平集函数正则性维持项、长度惩罚项、加权面积约束项，η,μ,ν,γ为权重系数；

①区域能量拟合项E_R(φ,f₁,f₂)定义为：

$E_R(\mathrm{\&phi;},f_1,f_2)=\underset{i=1}{\overset{2}{\&Sigma;}}{\mathrm{\&beta;}}_i\&Integral;(\&Integral;k_{\mathrm{\&sigma;}}(x-y\left){\left|I\left(y\right)-f_i\left(x\right)\right|}^2{M}_i^{\mathrm{\&epsiv;}}\right(\mathrm{\&phi;}\left(y\right)\left)dy\right)dx---\left(5\right)$

式(5)中k_σ(x-y)为高斯核函数，σ是高斯核尺寸参数，用以控制局部区域灰度值拟合范围；f_i(x),i＝1,2为局部区域灰度值拟合函数，通过标准梯度下降流法求解能量函数E(φ,f₁,f₂)得到：

$f_i\left(x\right)=\frac{k_{\mathrm{\&sigma;}}\left(x\right)*\&lsqb;M_i^{\mathrm{\&epsiv;}}\left(\mathrm{\&phi;}\right(x\left)\right)I\left(x\right)\&rsqb;}{k_{\mathrm{\&sigma;}}\left(x\right)*M_i^{\mathrm{\&epsiv;}}\left(\mathrm{\&phi;}\right(x\left)\right)},i=1,2---\left(6\right)$

式(6)β_i,i＝1,2为区域能量拟合项权重系数；φ为Lipschitz水平集函数，为区域函数，表示曲线内部和外部，定义H^ε(φ)为具有正则性的Heaviside函数，其导数为Dirac函数δ^ε(φ)；

$H^{\mathrm{\&epsiv;}}\left(\mathrm{\&phi;}\right)=\frac{1}{2}(1+\frac{2}{\mathrm{\&pi;}}arctan(\frac{\mathrm{\&phi;}}{\mathrm{\&epsiv;}}\left)\right)---\left(7\right)$

${\mathrm{\&delta;}}^{\mathrm{\&epsiv;}}\left(\mathrm{\&phi;}\right)=\frac{1}{\mathrm{\&pi;}}\frac{\mathrm{\&epsiv;}}{{\mathrm{\&epsiv;}}^2+{\mathrm{\&phi;}}^2}---\left(8\right)$

式(7)、式(8)中ε为尺寸参数；

由于作为水平集函数的初始化信息，即

②形状约束项E_S(φ)能量表达式定义为：

相应的水平集函数

式(9)中表示视网膜血管先验形状，即在t^*时刻

式(11)中d(x,C)表示视网膜图像中点x到轮廓曲线C的Euclidean距离，x∈C_outside表示点在轮廓线C外面，x∈C_inside表示点在轮廓线C内部，x∈C表示点在轮廓线C上；

③水平集函数正则性维持项P_R(φ)定义为：

$P_R\left(\mathrm{\&phi;}\right)=\frac{1}{2}{\&Integral;}_{\mathrm{\&Omega;}}{\left(\right|\&dtri;\mathrm{\&phi;}|-1)}^{2}dx---\left(12\right)$

④长度惩罚项L_P(φ)定义为：

$L_P\left(\mathrm{\&phi;}\right)={\&Integral;}_{\mathrm{\&Omega;}}|\&dtri;H^{\mathrm{\&epsiv;}}\left(\mathrm{\&phi;}\right)|dx---\left(13\right)$

⑤加权面积约束项A_R(φ)定义为：

A_R(φ)＝∫_Ωg(1-H^ε(φ))dx (14)

式(12)、式(13)和式(14)中Ω为图像全域，为非负单调递减函数，G为高斯函数，I_Rea等于I_B点乘I_Exten；在非血管的背景区域g＝1，面积约束项具有一个较大的惩罚，而在血管附近g→0，面积约束项只获得一个较小的惩罚；

(b)移除病灶和伪影，利用连通域面积Area和宽W、高H信息构建相应算子去除伪影与病灶；

当W/H和W×H同时满足：0.33<W/H<3.0和W×H<4Area时，则可判定整个连通域为病灶而非血管；当连通域面积Area小于50时，可判定整个连通域为伪影而非血管，从而得以提高视网膜血管的分割精度。

本发明在HRF和STARE、DRIVE数据库上的实验验证，准确率分别达到96.182％、95.034％和95.357％，优于现有分割方法。同时克服了在相邻血管处、血管交叉处和微血管处其它方法的不足，使分割出的血管结构最为接近金标准和血管真实尺寸值。

附图说明

图1为本发明实施例中的视网膜图像的绿色通道图像I。

图2为本发明实施例中自动获得视网膜的“掩模”。

图3为本发明实施例中消除了血管中心亮线的视网膜边缘扩展图像I_Exten。

图4为本发明实施例中视网膜图像背景估计I_Backg。

图5为本发明实施例中增强的视网膜血管图像I_enhance。

图6为本发明实施例中灰度值线性拉伸后获得的图像I_GrayS。

图7为本发明实施例中获得血管响应图像I_Resp。

图8为本发明实施例中选取阈值c₁获得的二值化图像

图9为本发明实施例中使用尺寸为2×2的模板对作开运算和删除连通域面积小于100的伪影，获得的二值化图像I_B。

图10为本发明实施例中的12_h视网膜图像伪影和病灶。

图11为本发明实施例中的12_h移除伪影和病灶后分割图像。

具体实施方式

下面结合具体实施方式，进一步阐述本发明。

实验说明：本发明的应用所涉及实施例数据来自于HRF数据库中第12个正常人(12_h)的视网膜图像。

本实施例包括三个步骤：视网膜血管图像预处理、血管图像粗略分割和血管图像精细分割。

具体描述如下：

1、视网膜血管图像预处理

(a)选取视网膜图像的绿色通道图像I，利用测地线活动轮廓(Geodesic activecontours，GAC)模型自动获得视网膜的“掩模”，如图2所示。

(b)利用上一步骤(a)获取的视网膜“掩模”信息，对图1做基于镜像对称的边缘扩展，边缘扩展的尺寸等于下一步骤(c)中的高斯模板尺寸。然后使用板尺寸为5×5形态学算子，对图像作开运算，消除血管中心亮线，可防止含有区域能量拟合项的水平集模型产生错误分割，从而获得消除了血管中心亮线的视网膜边缘扩展图像I_Exten，如图3所示。

(c)构建尺寸为100×100的高斯模板，对图像I_Exten做卷积运算，估计视网膜图像背景I_Backg，如图4所示，并与之做减法运算和与“掩模”做点乘运算，从而获得增强的视网膜血管图像I_enhance，如图5所示。然后把图像I_enhance的灰度值线性拉伸到[0-255]，获得图像I_GrayS，如图6所示。

2、视网膜血管图像粗略分割

视网膜血管图像预粗略分割主要包括以下两个步骤：

(a)构建血管响应函数

根据Hessian矩阵对图像不同点(x,y)的特征值描述，构建用以区别血管与背景的描述算子。本发明构建的视网膜血管描述算子J和血管响应函数ν_σ(x,y)分别定义为：

J＝(|λ₁|-|λ₂|)²(|λ₁|+|λ₂|)² (15)

$v_{\mathrm{\&sigma;}}(x,y)=\left\{\begin{array}{cc}0& {\mathrm{\&lambda;}}_2<0\\ (1-\exp (-\frac{R_{\mathrm{\&beta;}}^2}{2{\mathrm{\&beta;}}^2}\left)\right)(1-\exp (-\frac{J}{\sqrt{J_{\max }}}\left)\right)& {\mathrm{\&lambda;}}_2>=0\end{array}\right.---\left(16\right)$

式(15)、式(16)中1/R_β＝λ₂/λ₁，用于描述血管，λ₁和λ₂分别为二维Hessian矩阵的特征值。β为调整血管明暗的尺寸参数，这里取β＝0.5；σ是高斯核函数尺寸，J_max是J的最大值。上式(16)能最大化血管像素的响应，使血管响应函数ν_σ(x,y)→1_，而在非血管背景平坦区域上使血管响应函数ν_σ(x,y)→0，即ν_σ(x,y)∈[0,1]。

由于视网膜血管图像尺寸不一，需在多尺度下计算血管响应函数ν_σ(x,y)，最后取各尺度下的最大响应值，其定义如下：

$v(x,y)=\underset{{\mathrm{\&sigma;}}_{\min }\&le;\mathrm{\&sigma;}\&le;{\mathrm{\&sigma;}}_{max}}{{\mathrm{maxv}}_{\mathrm{\&sigma;}}(x,y)}---\left(17\right)$

式(17)中σ_min和σ_max为感兴趣的眼底视网膜血管的最小和最大尺寸。这里σ∈[1,11]，步长取2。

(b)二值化处理

利用式(16)和式(17)计算ν(x,y)，获得血管响应图像I_Resp(如图7)后进行二值化处理。本发明先选取较小的阈值c₁＝0.00001获得二值化图像如图8所示，并作为水平集函数φ的初始化。然后选取略大的阈值c₂＝0.0001获得二值化图像并使用尺寸为2×2的模板对作开运算和删除连通域面积小于100的伪影，获得二值化图像I_B，如图9所示。

3、视网膜血管图像精细分割

视网膜血管图像预精细分割主要包括以下两个步骤：

(a)构建视网膜图像血管分割的水平集模型

本发明利用粗略去除血管后视网膜背景方差设置形状约束项、面积约束项、长度惩罚项的权重系数，进一步克服图像噪声和固定时间开销。视网膜图像血管分割的水平集模型E(φ,f₁,f₂)定义如下：

E(φ,f₁,f₂)＝E_R(φ,f₁,f₂)+ηE_S(φ)+μP_R(φ)+νL_P(φ)+γA_R(φ) (18)

式(18)右边分别为局部区域能量拟合项、形状约束项、水平集函数正则性维持项、长度惩罚项、加权面积约束项。η,μ,ν,γ为权重系数。

①区域能量拟合项E_R(φ,f₁,f₂)定义为：

$E_R(\mathrm{\&phi;},f_1,f_2)=\underset{i=1}{\overset{2}{\&Sigma;}}{\mathrm{\&beta;}}_i\&Integral;(\&Integral;k_{\mathrm{\&sigma;}}(x-y\left){\left|I\left(y\right)-f_i\left(x\right)\right|}^2{M}_i^{\mathrm{\&epsiv;}}\right(\mathrm{\&phi;}\left(y\right)\left)dy\right)dx---\left(19\right)$

式(19)中k_σ(x-y)为高斯核函数，σ为高斯核尺寸参数，用以控制局部区域灰度值拟合范围；f_i(x),i＝1,2为局部区域灰度值拟合函数，通过标准梯度下降流法求解能量函数E(φ,f₁,f₂)得到

$f_i\left(x\right)=\frac{k_{\mathrm{\&sigma;}}\left(x\right)*\&lsqb;M_i^{\mathrm{\&epsiv;}}\left(\mathrm{\&phi;}\right(x\left)\right)I\left(x\right)\&rsqb;}{k_{\mathrm{\&sigma;}}\left(x\right)*M_i^{\mathrm{\&epsiv;}}\left(\mathrm{\&phi;}\right(x\left)\right)},i=1,2---\left(20\right)$

β_i,i＝1,2为区域能量拟合项权重系数，这里β₁＝β₂＝1；φ为Lipschitz水平集函数，为区域函数，表示曲线内部和外部。定义 H^ε(φ)为具有正则性的Heaviside函数，其导数为Dirac函数δ^ε(φ)。

$H^{\mathrm{\&epsiv;}}\left(\mathrm{\&phi;}\right)=\frac{1}{2}(1+\frac{2}{\mathrm{\&pi;}}arctan(\frac{\mathrm{\&phi;}}{\mathrm{\&epsiv;}}\left)\right)---\left(21\right)$

${\mathrm{\&delta;}}^{\mathrm{\&epsiv;}}\left(\mathrm{\&phi;}\right)=\frac{1}{\mathrm{\&pi;}}\frac{\mathrm{\&epsiv;}}{{\mathrm{\&epsiv;}}^2+{\mathrm{\&phi;}}^2}---\left(22\right)$

其中ε为尺寸参数，这里取ε＝1。

由于作为水平集函数的初始化信息，即

②形状约束项E_S(φ)能量表达式定义为：

相应的水平集函数

式(23)中表示视网膜血管先验形状，即在t^*时刻

其中d(x,C)表示视网膜图像中点x到轮廓曲线C的Euclidean距离，x∈C_outside表示点在轮廓线C外面，x∈C_inside表示点在轮廓线C内部，x∈C表示点在轮廓线C上。

③平集函数正则性维持项P_R(φ)定义为：

$P_R\left(\mathrm{\&phi;}\right)=\frac{1}{2}{\&Integral;}_{\mathrm{\&Omega;}}{\left(\right|\&dtri;\mathrm{\&phi;}|-1)}^{2}dx---\left(26\right)$

④长度惩罚项L_P(φ)定义为：

$L_P\left(\mathrm{\&phi;}\right)={\&Integral;}_{\mathrm{\&Omega;}}|\&dtri;H^{\mathrm{\&epsiv;}}\left(\mathrm{\&phi;}\right)|dx---\left(27\right)$

⑤加权面积约束项A_R(φ)定义为：

A_R(φ)＝∫_Ωg(1-H^ε(φ))dx (28)

式(26)、式(27)和式(28)中Ω为图像全域，为非负单调递减函数，G为高斯函数，I_Rea等于I_B点乘I_Exten。在非血管的背景区域g＝1，面积约束项具有一个较大的惩罚，而在血管附近g→0，面积约束项只获得一个较小的惩罚。

在视网膜血管分割实验中，时间步长t＝0.1，正则化项系数μ＝1、长度惩罚项系数ν＝ms、面积约束项系数γ＝nlog_as、形状约束项系数η＝ls²(s²表示粗略去除血管后视网膜背景方差)，本实施例中m、n、a、l分别取5、1、7、1。

(b)移除病灶和伪影

本发明利用连通域面积Area和宽、高(W,H)信息构建相应算子去除伪影与病灶。

本实施例构建8连通域模板标记水平集分割结果，从而获得每个连通域的面积Area和宽W、高H信息。

当W/H和W×H同时满足：

算子①：0.4<W/H<2.5

算子②：W×H<3.5Area

则整个连通域判定为病灶(非血管)。

算子③：Area<30，则整个连通域判定为伪影(非血管)。

利用算子①、算子②和算子③移除伪影和部分病灶，如图10所示，以提高视网膜血管的分割精度，如图11所示。

Claims (4)

1.一种融合形状先验的水平集视网膜血管图像分割方法，包括三个步骤：

步骤1，视网膜血管图像预处理：在对视网膜做边缘扩展基础上采用高斯卷积估计视网膜背景，以消除视网膜原始边界处的灰度突变，获得了一个背景更加均匀的血管增强图像；

步骤2，视网膜血管图像粗略分割：利用改进的血管响应函数和阈值得到视网膜血管的粗分割图像，并以此信息初始化水平集函数和构建形状约束项，以克服水平集模型对初始化敏感以及区域能量拟合项对噪声敏感的问题；

步骤3，视网膜血管图像精细分割：应用连通域面积Area和宽W、高H构建几何算子，消除连通域面积较小的伪影和病灶，获得高准确率的视网膜血管分割图像。

2.根据权利要求1所述的一种融合形状先验的水平集视网膜血管图像分割方法，其特征是：所述的步骤1包括以下三个子步骤：

(a)选取视网膜图像的绿色通道图像I，利用测地线活动轮廓模型自动获得视网膜的“掩模”；

(b)在获取的“掩模”的基础上对视网膜图像作边缘扩展，获得图像I_Exten；就是采用形态学算子，对图像作开运算，消除血管中心亮线，从而获得消除了血管中心亮线的视网膜边缘扩展图像I_Exten；

(c)构建高斯模板，所述的高斯模板尺寸为原始图像尺寸除以14后取整数，小数点后四舍五入；对图像I_Exten做卷积运算，估计视网膜图像背景I_Backg，并与I_Exten做减法运算和与“掩模”做点乘运算，从而获得增强的视网膜血管图像I_enhance，然后把图像I_enhance的灰度值线性拉伸到[0-255]，获得图像I_GrayS。

3.根据权利要求1所述的一种融合形状先验的水平集视网膜血管图像分割方法，其特征是：所述的步骤2包括以下两个子步骤：

(a)构建血管响应函数，根据Hessian矩阵对图像不同点(x,y)的特征值描述，构建用以区别血管与背景的描述算子；构建的视网膜血管描述算子J和血管响应函数ν_σ(x,y)分别定义为：

J＝(|λ₁|-|λ₂|)²(|λ₁|+|λ₂|)² (1)

${\mathrm{\&nu;}}_{\mathrm{\&sigma;}}(x,y)=\left\{\begin{array}{cc}0& {\mathrm{\&lambda;}}_2<0\\ (1-\exp (-\frac{R_{\mathrm{\&beta;}}^2}{2{\mathrm{\&beta;}}^2}\left)\right)(1-\exp (-\frac{J}{\sqrt{J_{\max }}}\left)\right)& {\mathrm{\&lambda;}}_2>=0\end{array}\right.---\left(2\right)$

式(1)、式(2)中1/R_β＝λ₂/λ₁，λ₁和λ₂分别为二维Hessian矩阵的特征值；β为调整血管明暗的尺寸参数，σ是高斯核函数尺寸，J_max是J的最大值；

由于视网膜血管图像尺寸不一，需在多尺度下计算血管响应函数ν_σ(x,y)，最后取各尺度下的最大响应值，其定义如下：

$\mathrm{\&nu;}(x,y)=\underset{{\mathrm{\&sigma;}}_{min}\&le;\mathrm{\&sigma;}\&le;{\mathrm{\&sigma;}}_{max}}{{\mathrm{max\&nu;}}_{\mathrm{\&sigma;}}(x,y)}---\left(3\right)$

式(3)中σ_min和σ_max为感兴趣的眼底视网膜血管的最小和最大尺寸；

(b)二值化处理，利用式(2)和式(3)计算ν(x,y)，获得血管响应图像I_Resp后进行二值化处理；先选取较小的阈值c₁获得二值化图像并作为水平集函数φ的初始化；然后选取略大的阈值c₂获得二值化图像并使用2×2的模板对作开运算和删除连通域面积小于100的伪影，获得二值化图像I_B去构建形状约束项。

4.根据权利要求1所述的一种融合形状先验的水平集视网膜血管图像分割方法，其特征是：所述的步骤3包括以下两个子步骤：

(a)构建视网膜图像血管分割的水平集模型，利用粗略去除血管后视网膜背景方差设置形状约束项、面积约束项、长度惩罚项的权重系数，进一步克服图像噪声和固定时间开销；视网膜图像血管分割的水平集模型E(φ,f₁,f₂)定义如下：

E(φ,f₁,f₂)＝E_R(φ,f₁,f₂)+ηE_S(φ)+μP_R(φ)+νL_P(φ)+γA_R(φ) (4)

式(4)右边分别为局部区域能量拟合项、形状约束项、水平集函数正则性维持项、长度惩罚项、加权面积约束项，η,μ,ν,γ为权重系数；

①区域能量拟合项E_R(φ,f₁,f₂)定义为：

$E_R(\mathrm{\&phi;},f_1,f_2)=\underset{i=1}{\overset{2}{\mathrm{\&Sigma;}}}{\mathrm{\&beta;}}_i\&Integral;(\&Integral;k_{\mathrm{\&sigma;}}(x-y\left){\left|I\left(y\right)-f_i\left(x\right)\right|}^2{M}_i^{\mathrm{\&epsiv;}}\right(\mathrm{\&phi;}\left(y\right)\left)dy\right)dx---\left(5\right)$

式(5)中k_σ(x-y)为高斯核函数，σ是高斯核尺寸参数，用以控制局部区域灰度值拟合范围；f_i(x),i＝1,2为局部区域灰度值拟合函数，通过标准梯度下降流法求解能量函数E(φ,f₁,f₂)得到：

$f_i\left(x\right)=\frac{k_{\mathrm{\&sigma;}}\left(x\right)*\&lsqb;M_i^{\mathrm{\&epsiv;}}\left(\mathrm{\&phi;}\right(x\left)\right)I\left(x\right)\&rsqb;}{k_{\mathrm{\&sigma;}}\left(x\right)*M_i^{\mathrm{\&epsiv;}}\left(\mathrm{\&phi;}\right(x\left)\right)},i=1,2---\left(6\right)$

式(6)β_i,i＝1,2为区域能量拟合项权重系数；φ为Lipschitz水平集函数，为区域函数，表示曲线内部和外部，定义H^ε(φ)为具有正则性的Heaviside函数，其导数为Dirac函数δ^ε(φ)；

$H^{\mathrm{\&epsiv;}}\left(\mathrm{\&phi;}\right)=\frac{1}{2}(1+\frac{2}{\mathrm{\&pi;}}\arctan (\frac{\mathrm{\&phi;}}{\mathrm{\&epsiv;}}\left)\right)---\left(7\right)$

${\mathrm{\&delta;}}^{\mathrm{\&epsiv;}}\left(\mathrm{\&phi;}\right)=\frac{1}{\mathrm{\&pi;}}\frac{\mathrm{\&epsiv;}}{{\mathrm{\&epsiv;}}^2+{\mathrm{\&phi;}}^2}---\left(8\right)$

式(7)、式(8)中ε为尺寸参数；

由于作为水平集函数的初始化信息，即

②形状约束项E_S(φ)能量表达式定义为：

相应的水平集函数

式(9)中表示视网膜血管先验形状，即在t^*时刻

式(11)中d(x,C)表示视网膜图像中点x到轮廓曲线C的Euclidean距离，x∈C_outside表示点在轮廓线C外面，x∈C_inside表示点在轮廓线C内部，x∈C表示点在轮廓线C上；

③水平集函数正则性维持项P_R(φ)定义为：

$P_R\left(\mathrm{\&phi;}\right)=\frac{1}{2}{\&Integral;}_{\mathrm{\&Omega;}}{\left(\right|\&dtri;\mathrm{\&phi;}|-1)}^{2}dx---\left(12\right)$

④长度惩罚项L_P(φ)定义为：

$L_P\left(\mathrm{\&phi;}\right)={\&Integral;}_{\mathrm{\&Omega;}}|\&dtri;H^{\mathrm{\&epsiv;}}\left(\mathrm{\&phi;}\right)|dx---\left(13\right)$

⑤加权面积约束项A_R(φ)定义为：

A_R(φ)＝∫_Ωg(1-H^ε(φ))dx (14)

式(12)、式(13)和式(14)中Ω为图像全域，为非负单调递减函数，G为高斯函数，I_Rea等于I_B点乘I_Exten；在非血管的背景区域g＝1，面积约束项具有一个较大的惩罚，而在血管附近g→0，面积约束项只获得一个较小的惩罚；

(b)移除病灶和伪影，利用连通域面积Area和宽W、高H信息构建相应算子去除伪影与病灶；

当W/H和W×H同时满足：0.33<W/H<3.0和W×H<4Area时，则可判定整个连通域为病灶而非血管；当连通域面积Area小于50时，可判定整个连通域为伪影而非血管，从而得以提高视网膜血管的分割精度。