CN104899862A - 基于全局与局部阈值的视网膜血管分割算法 - Google Patents

Abstract

本发明公开了一种基于全局与局部阈值的视网膜血管分割算法。该方法包括以下步骤：首先，提取视网膜图像中目标与背景对比度较高的绿色通道图像并采用多尺度多方向滤波方法来增强图像；其次，在增强图像上分别使用全局阈值二维最大熵与局部阈值移动平均算法对图像进行分割以获取视网膜的主血管与细小血管部分；最后，通过区域连通性的判断，将视网膜的主血管与细小血管进行结合，分割出最终的血管网络。本算法分割效果好，考虑到使用单一阈值方法无法有效提取整体血管网络，有效结合了全局阈值二维最大熵与局部阈值移动平均的优点，分割得到的血管在细小部分更为丰富，具有较好的连通性。本发明适用于正常及病变的视网膜图像的血管分割。

Description

基于全局与局部阈值的视网膜血管分割算法

技术领域

本发明涉及图像分割技术，尤其涉及一种基于全局与局部阈值的视网膜血管分割算法。

背景技术

综合分析现有的视网膜血管分割算法，较为简单常用的方法是阈值分割法，其中阈值分割又可分为全局阈值与局部阈值两种。基于全局阈值的方法中，其阈值的选取多依靠于灰度直方图，常用的方法有最大类间方差法(OTSU)和最大熵法等。针对视网膜血管分割方面，姚畅在论文《一种基于改进的PCNN的视网膜血管树提取方法》提出了一种基于改进的PCNN的视网膜血管树提取方法来提取视网膜血管网络，它根据视网膜血管网络的结构特征，对增强后的眼底图像运用IPCNN的动态点火特性分割出增强图像的血管网络，以提取出最终的血管树。该方法对整幅图像使用固定的全局阈值，如果图像中血管与背景的对比度十分明显，此分割方法相当有效。但是由于血管越往末梢越细，与背景的对比度也越来越弱，使用全局阈值只能得到血管的主干部分，血管的细小末梢部分难以分割。基于局部阈值的分割方法是对图像中的不同区域采用不同的阈值，其阈值的选取一般基于图像的局部统计信息，如局部方差、局部对比度以及曲 面拟合阈值等，又称为自适应阈值。针对视网膜血管分割方面，黄琳等在论文《视网膜图像中的血管自适应提取》提出了一种视网膜图像中的血管自适应提取的方法来提取视网膜血管网络，它通过视网膜图像进行分区然后计算每一区域满足梯度要求的像素点个数求取合适的局部阈值，最终实现血管的分割。该方法可以保证计算得到平均误差最小意义下的最优阈值，但是由于视网膜图像中背景灰度并不是非常均匀，所以阈值化得到的二值血管图像中包含很多碎片。

发明内容

本发明要解决的技术问题在于针对现有技术中的缺陷，提供一种基于全局与局部阈值的视网膜血管分割算法,本方法通过对增强后的视网膜图像分别进行全局与局部阈值分割，然后通过区域连通性的判断将二者结果进行融合，消除噪声，获取最终的血管分割结果。

本发明解决其技术问题所采用的技术方案是：一种基于全局与局部阈值的视网膜血管分割算法，包括以下步骤：

S1:对视网膜图像进行图像增强；

S2:对增强后的图像采用二维最大熵阈值法确定全局阈值，然后分割获得主血管部分的图像；

S3:对增强后的图像采用移动平均方法确定局部阈值，然后分割获得细小 血管部分的图像；

S4:对主血管部分和细小血管部分的图像进行融合，分割出最终的血管网络。

按上述方案，所述步骤S1)中的增强方法如下：

S11:对输入的原始RGB视网膜图像进行预处理，分解成红、绿、蓝三通道图像，提取视网膜图像的绿色通道图像；

S12:对绿色通道图像采用多尺度多方向滤波方法进行图像增强。

按上述方案，所述步骤S12)中的多尺度多方向滤波方法如下：

S121:输入绿通道图像，产生像素矩阵T；对于绿通道图像的每一个像素T_ab，首先初始化空间尺度σ的范围[a₁,a₂]、迭代步长Step、滤波器的增强因子β和c；

如果σ没有满足停止条件，那么根据血管相似滤波器计算增强滤波的输出值Z_o，接着迭代σ，判断σ是否满足停止条件；

如果σ满足停止条件，则选取最大增强滤波输出值作为该像素的增强因子；

S122:对于绿通道图像的每一个像素T_ab，首先初始化方向因子θ的范围[a₁,a₂]、迭代步长Step；

如果θ没有满足停止条件，求得相应的高斯核模板并与T进行卷积，得卷积结果Z_θ；迭代θ的值，判断θ是否满足停止条件；

如果满足停止条件，则选取所有方向相应像素的最大值作为该像素的增强 因子；

S123:对比分析每一个像素f_ij的相应位置多尺度增强结果f₁和多方向增强结果f₂的值，选取f₁和f₂中最大值作为f_ij的结果，从而得到最终的增强结果。

按上述方案，所述步骤S3)中局部阈值的确定方法具体为：

移动平均的扫描方式为Z字形逐行扫描；

令z_k+1表示步骤k+1中扫描序列中遇到的点的灰度，则z_k+1处的移动平均灰度为其中k≥n-1；n表示用于计算平均的点数，且m(1)＝z₁/n；

当k<n-1时，移动平均灰度

对于图像中的每个点都计算移动平均，则可变局部阈值的通用形式为：T_xy＝bm_xy；其中，b是常数，m_xy是在输入图像中的点(x,y)处计算得到的移动平均。

按上述方案，所述步骤S4)中的融合方法具体如下：将二维最大熵结果中的主血管作为整个血管网络的主干，通过对区域连通性进行判断，即当移动平均中像素为细小血管图像中血管像素，并且在大小为a×a的邻域窗口中有一个像素在相应的主血管图像中为血管像素时，即认为该点像素为血管，然后循环遍历整个移动平均图像，直至判定完成。

本发明产生的有益效果是：

1.本发明解决了使用单一阈值方法无法有效提取整体血管网络的问题，利用全局阈值二维最大熵提取视网膜图像中主血管网络，同时利用局部阈值移动平均提取视网膜图像中细小血管部分，然后通过区域连通性的判断将以上两个分割结果进行融合，从而实现血管的最终分割。

2.本发明分割效果好，在保证主血管完整性的同时不丢失与主血管相连的细小血管。

附图说明

下面将结合附图及实施例对本发明作进一步说明，附图中：

图1为本发明实施例的方法流程图；

图2为本发明实施例的预处理结果图；

图3为本发明实施例的血管分割结果。

图4为眼底照相机拍摄得到的一幅正常视网膜图像；

图5为图4视网膜图像使用本发明方法的分割结果；

图6为眼底照相机拍摄得到的一幅病变视网膜图像。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例， 对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。

图1为本发明基于全局与局部阈值的视网膜血管分割算法的流程图，如图1所示，本方法包括以下步骤：

步骤S1:对视网膜图像进行图像增强；

对输入的原始RGB视网膜图像进行预处理，分解成红、绿、蓝三通道图像，由于绿通道的图像中血管与背景的对比度最高，因此在绿通道图像上采用多尺度多方向滤波方法来增强图像。

输入绿通道图像，产生像素矩阵T。对于T中的每一个像素T_ab，首先初始化空间尺度σ的范围[a₁,a₂]、迭代步长Step、滤波器的增强因子β和c。

如果σ没有满足停止条件，那么根据公式计算元素T_ab与高斯函数二阶微分的卷积。然后，根据血管相似滤波器 计算增强滤波的输出值Z_o，接着迭代σ，判断σ是否满足停止条件(σ≤a₂)。

如果σ满足停止条件，则根据公式选取最大增强滤波输出值Z_max作为该像素的尺度增强因子。

在本发明实施例中，设定空间尺度σ的范围为[0.2,3]，迭代步长Step为0.1。

输入绿通道图像，产生像素矩阵T。对于T中的每一个像素T_ab，首先初始化方向因子θ的范围[a₁,a₂]、迭代步长Step。如果θ没有满足停止条件(θ≤a₂)，那么对于每个θ，根据公式c′_i(x,y)＝c_i(x,y)-a_i求得相应的高斯核模板。其中，c_i(x,y)为相应高斯核矩阵系数，可由公式求得，u为旋转后相应的纵坐标，a_i为相应的高斯核矩阵的平均值。

将计算得到的高斯核模板与T进行卷积，得卷积结果Z_θ。迭代θ的值，判断θ是否满足停止条件。如果满足停止条件，对于T中的每一个像素T_ab根据公式 取所有方向相应像素的最大值Z_max作为该像素的方向增强因子，其中g(i，j)为绿通道图像在(i，j)上的像素值。

在本发明实施例中，设定方向因子θ的范围为[0,165]，迭代步长Step为30。

对比分析每一个像素f_ij的相应位置多尺度增强结果f₁和多方向增强结果f₂的值，选取f₁和f₂中最大值作为f_ij的结果，即最终的增强结果可由公式 表示。

步骤S2:利用第一步得到的增强后的图像采用全局阈值中的二维最大熵阈值法分割主血管部分；

图像的灰度取值范围为{0,1,2,…,L-1}，假设分割阈值矢量为(S,T)，将图像分成目标A与背景C，它们的概率分布分别是： 其中， $p_{st}={\mathrm{\&Sigma;}}_{i=0}^s{\mathrm{\&Sigma;}}_{j=0}^t{p}_{ij},p_{ij}=g_{ij}/(M*N),$ g_ij为相应灰度与领域平均灰度对出现的频数，M*N为图像的大小。

目标A与背景C的熵分别被定义为： $H_C(s,t)=-{\mathrm{\&Sigma;}}_{i=s+1}^{L-1}{\mathrm{\&Sigma;}}_{j=t+1}^{L-1}\frac{p_{ij}}{1-p_{st}}\ln \frac{p_{ij}}{1-p_{st}}.$

为了获取图像中目标和背景的最大信息量，根据公式 求得当总熵H(s,t)最大时的灰度，对(S,T)即为所求的最佳阈值。

步骤S3:利用第一步得到的增强后的图像采用局部阈值中的移动平均方法分割细小血管部分；

移动平均的扫描方式是以Z字形模式逐行执行的。令z_k+1表示步骤k+1中扫描序列中遇到的点的灰度，这个新点处的移动平均(平均灰度)可由  $m(k+1)=\frac{1}{n}{\mathrm{\&Sigma;}}_{i=k+2-n}^{k+1}{z}_i=m\left(k\right)+\frac{1}{n}(z_{k+1}-z_{k-n})$ 计算求得。其中，n表示用于计算平均的点数，且m(1)＝z₁/n。第一个表达式在k≥n-1时有效，当k<n-1时，平均是使用可用点形成的。类似地，第二个表达式在k≥n+1时有效。

对于图像中的每个点都计算移动平均，可变局部阈值的通用形式为：T_xy＝bm_xy；其中，b是常数，m_xy是在输入图像中的点(x,y)处使用 $m(k+1)=$ $\frac{1}{n}{\mathrm{\&Sigma;}}_{i=k+2-n}^{k+1}{z}_i=m\left(k\right)+\frac{1}{n}(z_{k+1}-z_{k-n})$ 计算得到的移动平均。

在本发明的实施例中，计算平均的点数n为20，b为0.5。

步骤S4:对第二步和第三步得到的图像进行融合，分割出最终的血管网络。

基于二维最大熵分割出的血管是血管增强图像中的主血管部分，对于图像中的小血管部分不能进行很好地分割。基于移动平均提取出了大部分的细小血管和血管末梢部分，但同时也包含了部分背景与噪声。为了区分开移动平均中的血管与噪声部分，将二维最大熵结果中的主血管作为整个血管网络的主干，通过对区域连通性进行判断，即当移动平均中像素为血管像素，并且在大小为a×a的邻域窗口中有一个像素在相应的主血管图像中为血管像素时，即认为该点像素为血管，然后循环遍历整个移动平均图像，直至判定完成。这样，就将移动平均中位于主血管周围的小血管提取出来了，而移动平均中与主血管没有连通性的噪声部分则被清除。

在本发明的实施例中，邻域窗口大小为5×5。

为验证本发明的有效性和实用性，接下来在在国际公认的视网膜图像数据库Drive上进行了实验。这个数据库提供20幅视网膜图像用于算法测试，两个手动分割的血管网络数据集作为参考。选择数据库中的1幅图像进行处理，如图2、3所示。其中，图2(a)为原始彩色视网膜图像，图2(b)为绿色通道图像，图2(c)为增强后的视网膜图像，即视网膜图像的预处理结果；图3(a) 为二维最大熵分割结果，图3(b)为移动平均分割结果，图3(c)为该视网膜图像的血管结果。由图3所示的实验结果可知，本发明能够有效地对视网膜图像的整体血管网络进行提取分割。

图4为眼底照相机拍摄得到的一幅正常视网膜图像；图5为图4视网膜图像使用本发明方法的分割结果；图6为眼底照相机拍摄得到的一幅病变视网膜图像。由图4至图6所示的结果可知，本发明方法适用于正常及病变的视网膜图像的血管分割。

应当理解的是，对本领域普通技术人员来说，可以根据上述说明加以改进或变换，而所有这些改进和变换都应属于本发明所附权利要求的保护范围。

Claims (6)

1.一种基于全局与局部阈值的视网膜血管分割算法，其特征在于，包括以下步骤：

S1:对视网膜图像进行图像增强；

S2:对增强后的图像采用二维最大熵阈值法确定全局阈值，然后分割获得主血管部分的图像；

S3:对增强后的图像采用移动平均方法确定局部阈值，然后分割获得细小血管部分的图像；

S4:对主血管部分和细小血管部分的图像进行融合，分割出最终的血管网络。

2.根据权利要求1所述的视网膜血管分割算法，其特征在于，所述步骤S1)中的增强方法如下：

S11:对输入的原始RGB视网膜图像进行预处理，分解成红、绿、蓝三通道图像，提取视网膜图像的绿色通道图像；

S12:对绿色通道图像采用多尺度多方向滤波方法进行图像增强。

3.根据权利要求2所述的视网膜血管分割算法，其特征在于，所述步骤S12)中的多尺度多方向滤波方法如下：

S121:输入绿通道图像，产生像素矩阵T，对于绿通道图像的每一个像素，首先初始化空间尺度σ的范围、迭代步长、滤波器的增强因子β和c；

如果σ没有满足停止条件，那么根据血管相似滤波器计算增强滤波的输出值Z_o，接着迭代σ，判断σ是否满足停止条件；

如果σ满足停止条件，则选取最大增强滤波输出值作为该像素的增强因子；

S122:对于绿通道图像的每一个像素，首先初始化方向因子θ的范围和迭代步长；

如果θ没有满足停止条件，求得相应的高斯核模板并与T进行卷积，得卷积结果Z_θ；迭代θ的值，判断θ是否满足停止条件；

如果满足停止条件，则选取所有方向相应像素的最大值作为该像素的增强因子；

S123:对比分析每一个像素f_ij的相应位置多尺度增强结果f₁和多方向增强结果f₂的值，选取f₁和f₂中最大值作为f_ij的结果，从而得到最终的增强结果。

4.根据权利要求1所述的视网膜血管分割算法，其特征在于，所述步骤S3)中局部阈值的确定方法具体为：

移动平均的扫描方式为Z字形逐行扫描；

令z_k+1表示步骤k+1中扫描序列中遇到的点的灰度，则z_k+1处的移动平均灰度为其中k≥n-1；n表示用于计算平均的点数，且m(1)＝z₁/n；

当k<n-1时，移动平均灰度

对于图像中的每个点都计算移动平均，则可变局部阈值的通用形式为：T_xy＝bm_xy；其中，b是常数，m_xy是在输入图像中的点(x,y)处计算得到的移动平均。

5.根据权利要求1所述的视网膜血管分割算法，其特征在于，所述步骤S4)中的融合方法具体如下：将二维最大熵结果中的主血管作为整个血管网络的主干，通过对区域连通性进行判断，即当移动平均中像素为细小血管图像中血管像素，并且在大小为a×a的邻域窗口中有一个像素在相应的主血管图像中为血管像素时，即认为该点像素为血管，然后循环遍历整个移动平均图像，直至判定完成。

6.根据权利要求5所述的视网膜血管分割算法，其特征在于，所述a取值为5。