基于非荧光眼底图像的视网膜血管自动分割方法

摘要

本发明涉及一种基于非荧光眼底图像的视网膜血管自动分割方法，包括以下步骤：第一步、对图像进行预处理，增强血管特征，弱化背景噪声；第二步、利用多尺度线性算子提取视网膜图像的形态学特征，并结合二维Gabor小波变换分割出精细的血管框架；第三步、分割出精细的血管框架之后，采用路径形态学滤波，将最大路径长度与血管区域的几何特征相结合，构造临接图；第四步、通过区域连通性分析和滞后阈值技术进行二值化过程，辅助分割微小血管，完成分割。与现有技术相比，本发明血管定位精度较高，粘连现象少，有效分割眼底图像的视网膜血管网络。

说明书

技术领域

本发明涉及图像处理技术领域，尤其是涉及一种基于非荧光眼底图像的视网膜血管自动分割方法。

背景技术

1989年，Chaudhuri等人在发表的“Detection of blood vessels in retinal images using two-dimensional matched filters”中最早提出了二维匹配滤波方法，该方法假定血管为一些等宽的线段，血管宽度介于2-10个像素，且血管横截面的灰度分布可以通过高斯曲线来近似，基于此构造12个不同方向的高斯模板来对图像进行匹配滤波，将最大响应结果作为滤波结果输出。该方法当时针对视网膜图像中的血管取得了出色的增强效果，但不足的是，当图像中摻杂与血管类似高斯分布的噪声时，噪声经过匹配滤波器后也会进行增强，与血管难以区分。此外，该方法没有考虑如何对滤波后的图像进行阈值化处理，无法与专家手动分割的标准图像进行比较。

Zana等人在发表的“Segmentation of vessel-like patterns using mathematical morphology and curvature evaluation”中巧妙地将血管的剖面曲率评价方法加入形态学算子，以提高血管分割的精度。其核心思想有两点：(1)形态学算子考虑了血管的形态学特征，如连续，分段和线性；(2)剖面曲率评价方法考虑了血管结构的特质：其剖面曲率是线性连贯的。文中，作者首先对图像进行预处理，平滑背景部分以降低噪声的影响；其次，采用类高斯分布的线性结构元素对图像中血管进行增强，并动态跟踪血管的剖面曲率，以此将血管和背景区别开。但该算法对待处理图像本身的质量要求较高，若存在较大的高光或者暗斑等光照不均匀区域，血管分割的精度会受到较大影响；同时，形态学算法的关键在于选取恰当的结构元素(菱形、圆形、矩形等)，若选取不当，血管形态会失准。

近年来，国际上研究较多且分割效果较好的是监督分割类方法。这类方法的特点是从专家手动分割的血管图像中提取特定的血管特征向量，以此作为训练集来优化分类器模型；当分类器模型优化到一定程度后，就可取得较好的分割效果。Niemeijer等人利用窗口估计提取眼底图像的绿色分量灰度值、局部梯度最大值和局部最大特征值作为特征向量进行分割。监督类分割方法利用了视网膜血管标准图像信息，分割结果令人满意，但需要大量手动专家分割的标准图像，分割效率不高。

目前相关文献提出的方法普遍存在微小血管分割精度不高和易出现粘连现象的问题。

发明内容

本发明的目的就是为了克服上述现有技术存在的缺陷而提供一种基于非荧光眼底图像的视网膜血管自动分割方法，作为一种非监督类分割算法，该方法血管定位精度较高，粘连现象少，有效分割眼底图像的视网膜血管网络。

本发明的目的可以通过以下技术方案来实现：

一种基于非荧光眼底图像的视网膜血管自动分割方法，包括以下步骤：

第一步、对图像进行预处理，增强血管特征，弱化背景噪声；

第二步、利用多尺度线性算子提取视网膜图像的形态学特征，并结合二维Gabor小波变换分割出精细的血管框架；

第三步、分割出精细的血管框架之后，采用路径形态学滤波，将最大路径长度与血管区域的几何特征相结合，构造临接图；

第四步、通过区域连通性分析和滞后阈值技术进行二值化过程，辅助分割微小血管，完成分割。

所述的第一步具体包括以下子步骤：

(1)图像进行预处理的过程包括对比度增强和视网膜边界生长；

(2)明采用CLAHE算法来改进G通道图像的局部对比度，期望呈现更多图像细节；当图像的直方图幅度超过设定阈值，CLAHE算法会对其进行裁剪限幅，将裁剪掉部分均匀分布到原直方图的其他部分；

(3)视网膜边界生长基于迭代式形态学膨胀操作。

所述的步骤(2)中的G通道图像划分为64个局部区域，分别计算直方图、累计分布函数以及对应的灰度变换函数，直方图幅度阈值设置为0.01。

所述的基于迭代式形态学膨胀操作具体过程为：

(a)通过形态学腐蚀操作来定位眼底图像的视网膜边界区域；

(b)迭代地采用形态学膨胀操作对视网膜边界区域向外扩充，在扩充过程中，对边界上的像素点进行均值滤波。

所述的第二步具体包括以下内容：

将线性算子扩展到多尺度，且该线性算子可以旋转以匹配不同方向的血管，不仅排除了中央反射的干扰，而且减少了血管分叉点或平行紧挨血管的分割误差；

假设输入的眼底图像数学定义为I(i，j)，L_s,k(i,j)为多尺度线性算子与输图像的卷积响应，对于每一个尺度，该线性算子都会旋转以匹配不同方向的血管，真正匹配时所得到的响应输出h_s(i,j)为：

$h_s(i,j)=\underset{k}{\max }{L}_{s,k}(i,j)-n(i,j)$

式中，尺度取值范围为s＝1,3,...,15，方向个数为k＝0,1,...,11，n(i,j)表示灰度响应的均值输出，基于窗口大小为w＝15计算的；

将大尺度和小尺度的优势结合起来，得到多尺度的线性算子响应H(i,j)，表达式如下：

$H(i,j)=\frac{\underset{s}{\Sigma }{h}_s(i,j)+I(i,j)}{n_s+1}$

式中，n_s为线性算子不同尺度的数量；

在H(i,j)的表达式中引入了局部区域方差参数α，，虽然实际输出响应增大，但α值也会增大，从而区分开背景与血管，H(i,j)的表达式更新为：

$H(i,j)=\left(\begin{array}{c}0,if(\alpha >t_{\alpha })\\ \frac{\underset{s}{\Sigma }{h}_s(i,j)+I(i,j)}{n_s+1}\end{array}\right)$

所述的第三步具体包括以下内容：

将形态学滤波方法分别应用于多尺度线性算子和二维Gabor小波变换后的结果图像，并选取设定的长度阈值和容忍度，有效滤除了噪声，改进了不连续血管的分割效果，为后面二值化过程奠定了良好基础。

所述的第四步具体包括以下内容：

采用滞后阈值方法对路径形态学滤波后图像进行二值化处理，具体为，设定低阈值t_l和高阈值t_h，来衡量图像中像素点对于血管边缘响应的幅值；若幅值超过高阈值t_h，将这些像素当作是血管区域；若幅值低于高阈值t_h，则分两种情况区别对待：(1)幅值低于低阈值t_l，不考虑；(2)幅值超过低阈值t_l，如果该像素点与周围的像素点有邻接关系，且周围的像素点已经被保留血管，则将这些像素点也看作是血管；不断迭代上述操作，直至所有的像素点都被标记为噪声或血管。

与现有技术相比，本发明具有以下优点：

1)针对现有血管分割方法的不足，尤其是血管中央反射和病变区域的分割问题，本文将线性算子扩展到多尺度，引入局部区域方差参数，能有效分割眼底图像的视网膜血管网络，且该线性算子可以旋转以匹配不同方向的血管，保留了较为完整的血管信息，受噪声干扰小；

2)使用二维Gabor小波变换和路径形态学滤波。利用Gabor小波提取视网膜图像的空间域特征、频域特征和形态学特征，能够准确地分析信号的方向特征，适合于提高眼底图像微小血管与背景的对比度。并采用路径形态学滤波构造临接图，能够灵活地适应血管略微弯曲和曲率较大的形态结构，较好地抑制了光照不均等噪声干扰

附图说明

图1是本发明方法的流程图。

图2是眼底图像预处理过程示意图。(a)原始图像；(b)G通道图像；(c)使用CLAHE算法改进的G通道图像；(d)CLAHE算法的迭代效果；(e)视网膜边界生长。

图3是视网膜边界增长过程示意图。

图4是多尺度线性算子分割结果。

图5是多尺度线性算子和二维Gabor小波变换对于细小血管提取效果的比较。(a)眼底图像局部区域放大部分；(b)多尺度线性算子输出响应；(c)二维Gabor小波变换在尺度因子值为1.5时的输出响应；(d)二维Gabor小波变换在尺度因子值为3时的输出响应。

图6是形态学滤波实验结果，左边为原图，右边为滤波后的结果。

图7是二值化过程。(a)滞后阈值；(b)长度滤波。

具体实施方式

下面结合附图和具体实施例对本发明进行详细说明。本实施例以本发明技术方案为前提进行实施，给出了详细的实施方式和具体的操作过程，但本发明的保护范围不限于下述的实施例。

实施例

本实施例包括以下步骤：

第一步、对图像进行预处理，增强血管特征，弱化背景噪声。

所述的预处理包括对比度增强和视网膜边界生长。

对比度增强主要基于限制对比度自适应直方图均衡化算法。R通道过度曝光，对比度低；B通道亮度较低，血管不易辨认；相比R和B通道，G通道图像血管与背景之间对比度最高，噪声较少。因此，我们选用G通道图像来进行后续处理。本发明采用CLAHE算法来改进G通道图像的局部对比度，期望呈现更多图像细节。和普通的自适应直方图均衡方法相比，CLAHE的特点在于其对比度限幅，从而克服了普通方法中存在的噪声过度放大问题。当图像的直方图幅度超过某个阈值，CLAHE算法会对其进行裁剪限幅，将裁剪掉部分均匀分布到原直方图的其他部分。重分布过程可能会导致图像的直方图幅度仍然超过阈值，如果这不是所希望的，可以迭代该操作直至满足要求。本发明中，我们将G通道图像划分为64个局部区域，分别计算直方图、累计分布函数以及对应的灰度变换函数，直方图幅度阈值设置为0.01

视网膜边界生长主要基于迭代式形态学膨胀操作。具体过程为：

(1)通过形态学腐蚀操作来定位眼底图像的视网膜边界区域。

(2)迭代地采用形态学膨胀操作对视网膜边界区域向外扩充，在扩充过程中，对边界上的像素点进行均值滤波。

第二步、利用多尺度线性算子提取视网膜图像的形态学特征，并结合二维Gabor小波变换分割出精细的血管框架。

线性算子在视网膜血管分割方面其主要原理在于，血管的灰度明显比背景低，通过计算均值和线性算子的最大响应，并将二者比较作差，即可区分出血管与背景。该方法的特点在于简单高效且分割出的血管结构轮廓清晰，而且限制少，适用于大多数应用场景。本发明将线性算子扩展到多尺度，且该线性算子可以旋转以匹配不同方向的血管，不仅排除了中央反射的干扰，而且减少了血管分叉点或平行紧挨血管的分割误差。

假设输入的眼底图像数学定义为I(i，j)，L_s,k(i,j)为多尺度线性算子与输入图像的卷积响应，对于每一个尺度，该线性算子都会旋转以匹配不同方向的血管，真正匹配时所得到的响应输出h_s(i,j)为：

$h_s(i,j)=\underset{k}{\max }{L}_{s,k}(i,j)-n(i,j)$

式中，尺度取值范围为s＝1,3,...,15，方向个数为k＝0,1,...,11，n(i,j)表示灰度响应的均值输出，基于窗口大小为w＝15计算的。小尺度线性算子在正确分割血管结构的同时引入了大量噪声，大尺度则有效抑制了噪声；但大尺度线性算子在一些复杂情况(如血管分叉点或血管之间挨得很近)下的分割结果不尽如人意，易出现粘连情况，小尺度则合理地规避了这一点。

本发明巧妙地将大尺度和小尺度的优势结合起来，得到多尺度的线性算子响应H(i,j)，表达式如下：

$H(i,j)=\frac{\underset{s}{\Sigma }{h}_s(i,j)+I(i,j)}{n_s+1}$

式中，n_s为线性算子不同尺度的数量。

我们观察到，在血管之间非常靠近的情况下，经典的线性算子方法分割结果会出现血管融合的现象。但血管之间靠的很近，只有一个像素宽的距离，且目标像素本身属于背景。由于线性算子方法在计算过程中会涉及到目标像素的领域，可能会导致实际的输出响应增大，有误判的风险。为了避免这种情况的出现，我们在H(i,j)的表达式中引入了局部区域方差参数α，上述情况下，虽然实际输出响应增大，但α值也会增大，从而区分开背景与血管。H(i,j)的表达式更新为：

$H(i,j)=\left(\begin{array}{c}0,if(\alpha >t_{\alpha })\\ \frac{\underset{s}{\Sigma }{h}_s(i,j)+I(i,j)}{n_s+1}\end{array}\right)$

多尺度线性算子获得了良好的血管特征检测效果，粘连现象少，有效滤除了光照不均等噪声。不足的是，由于部分细小血管对比度较低，其分割结果缺失了部分细节信息。这个问题用二维Gabor小波变换方法解决。

二维Gabor小波是方向性小波，能够准确地分析信号的向特征，适合于提高眼底图像微小血管与背景的对比度。因此，本文采用Gabor小波来增强OAD算子对微小血管结构的敏感度。本发明中的数学公式和符号参照Arneodo等人的论文”A wavelet-based method for multifractal image analysis.I.Methodology and test applications on isotropic and anisotropic random rough surfaces”。

$T_{\psi }(b,\theta ,a)={C_{\psi }}^{-1/2}\frac{1}{a}\int {\psi }^{*}(a^{-1}r-\theta (x-b\left)\right)I\left(x\right)d^{2}x$

ψ(x)＝exp(jk₀x)exp(-xAx/2)

r_θ(x)＝(xcosθ-ysinθ,xsinθ+ycosθ)

式中*表示共轭，C_ψ表示归一化常量，ψ表示要分析的小波，b为位移矢量，θ为旋转角度，a为尺度系数，k₀为复指数频矢量。A＝diag[η^-1/2,1]为一个对角矩阵，其中η取值较为关键，它控制在任何指定方向上的滤波器伸长。本发明参照Soares等人在“Retinal>0＝[0,2,5]。这样取值在增强细小血管特征的同时抑制了非血管噪声的放大。

小波变换系数T_ψ(b,θ,a)取值与实际小波的频率呈正相关，当眼底图像中血管的灰度分布特征匹配母小波的波形曲线时，T_ψ(b,θ,a)的模值达到峰值。其中，血管的脊线(ridge)定义为T_ψ(b,θ,a)取峰值时的位置坐标集合。

ridge(b,a)＝max_θ|T_ψ(b,θ,a_i)|

式中，i为尺度系数a的数量，a_i为各尺度系数的值，θ代表着母小波的旋转角。因为求得的ridge与尺度、方向有关，所以采用多尺度、多方向更加合理，即在图像的每个像素点选择多个值(本发明选择4个尺度，18个方向)，将小波变换系数，模最大值ridge(b,a)作为当前点的特征函数值输出，这样就得到了视网膜血管的特征函数图。

第三步、采用路径形态学滤波，将最大路径长度与血管区域的几何特征相结合，构造临接图。

路径形态学是在数学形态学的基础上发展而来，其核心思想是构造一组特定形态的结构元素，即路径，去匹配图像中对应的目标形状。其中，路径算子的构造非常关键，通常在邻接图上定义，其中邻接关系主要有四种：两种方向用于模拟水平和垂直，两种用于模拟对角线方向。路径算子主要有两个参数：路径长度和容忍度。路径长度决定了最终保留对象的大小，常以血管的分支长度作为阈值。容忍度表示在路径上允许目标对象像素不连续的像素个数。

通过设定路径算子长度阈值N_min和容忍度K，长度大于或等于N_min的区域会被保留下来，同时在路径上的任意位置，最多有K个像素缺失。

在彩色眼底图像中，血管局部表现为灰度均匀，且曲率变化平缓，呈狭长管状区域；同时由于眼底摄像技术复杂，眼底本身也会产生反射光干扰，导致实际获取的图像存在大面积的噪声，这可能会造成原本完整连续的血管中间部分有缺失，从而被截断为若干不连续的细小血管。本发明将形态学滤波方法分别应用于多尺度线性算子和二维Gabor小波变换后的结果图像，并选取合适的长度阈值和容忍度，有效滤除了噪声，较大地改进了不连续血管的分割效果，为后面二值化过程奠定了良好基础。

第四步、通过区域连通性分析和滞后阈值技术进行二值化过程，辅助分割微小血管。

本文采用滞后阈值方法对路径形态学滤波后图像进行二值化处理。其基本思路是，设定低阈值t_l和高阈值t_h，来衡量图像中像素点对于血管边缘响应的幅值。若幅值超过高阈值t_h，这些像素将毫无疑问被当作是血管区域；若幅值低于高阈值t_h，则分两种情况区别对待：(1)幅值低于低阈值t_l，不考虑；(2)幅值超过低阈值t_l，如果该像素点与周围的像素点有邻接关系(4-邻接或8-邻接)，且周围的像素点已经被保留血管，则将这些像素点也看作是血管。不断迭代上述操作，直至所有的像素点都被标记为噪声或血管。实验中，经反复验证发现，通过自适应Ostu算法选取t_h，取t_l为t_h的1/2时，二值化结果最佳。

由于视网膜血管在眼底图像上呈现连通的树状网络结构，可以通过区域连通性分析方法去除所有长度小于40个像素的分支，从而排除大部分残留非血管像素以提取更为完整、精确的血管骨架结构，得到最终的血管树结构

实施效果

依据上述步骤，对于我们收集的眼底图像进行分析，我们采用国际上公开的视网膜眼底图像数据库DRIVE和STARE中的图像进行实验比较进行实验。所有试验均在PC计算机上实现，该PC计算机的主要参数为：中央处理器Intel(R)Core(TM)2 Duo CPU E7500@2.93GHz，内存2GB。

结果显示，本发明对两个数据库(DRIVE和STARE)中图片进行视网膜血管自动分割的准确度分别为80.74％和68.87％。准确度方面明显优于其他算法。有效消除了视盘的影响，粘连现象少，同时能够保留更多的血管末梢细节。但发明方法在分割过程中也出现了部分图像细节的缺失，如血管形态方面，和专家手动分割图像对比，本发明算法在毛细血管的分割精度方面有所欠缺，存在过分割效应。原因在于，本文采用高斯核函数进行匹配滤波，虽然提高了血管对比度，但也在一定程度上模糊了血管边缘，对细小血管的定位不够准确。此外，分割结果中仍存在个别不连通现象，“类血管状”病变也不能够完全去除，有所改进。总体来说，本发明算法分割效果较好。

以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到各种等效的修改或替换，这些修改或替换都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以权利要求的保护范围为准。