视网膜血管图像的分割方法与流程

本发明涉及数字图像处理技术领域，尤其涉及一种视网膜血管图像的分割方法。

背景技术：

视网膜血管是人体全身血管中唯一可以无创直接观察到的血管，它的形状、管径、尺度、分支角度是否有变化，以及是否有增生、渗出，均可反映全身血管的病变。糖尿病、高血压、脑血管硬化等疾病均可导致视网膜血管发生一定的变化，因此视网膜血管的检测和分析对这些疾病的诊断和治疗在临床上具有重要的指导意义。然而，由于视网膜血管图像的灰度分布不均匀，目标血管与图像背景的对比度低，再加上图像噪声的污染，使视网膜血管的自动分割非常困难。

目前，现有技术中的一种视网膜血管图像的分割方法为：视网膜血管追踪方法。Tolias(1998)等应用模糊C均值方法，从视网膜圆盘处出发建立种子点，利用血管断面的一维模糊模型进行视网膜血管的追踪，通过判定所有种子点是否属于血管进行最终分割提取。这种方法能比较全面地描述视网膜血管网络结构。

上述视网膜血管追踪方法的缺点为：运算量大，而且对于血管的分支点以及对比度较低的血管分割不够精确。这种方法对种子点的选取比较敏感，算法往往在血管分支点处终止，丢失大量视网膜小血管。

现有技术中的另一种视网膜血管图像的分割方法为：分类器识别方法。该方法主要是通过对视网膜血管进行预处理操作，根据血管的不同特性建立特征空间，选择合适的分类器对样本进行训练，然后带入测试集中判定最后的结果。Staal(2004)采用了脊线检测方法，利用血管中心线邻域的特征向量进行有监督类识别。Soares(2006)首先抽取视网膜血管图像的多尺度2D Gabor小波特征，然后采用贝叶斯分类器对视网膜血管进行识别。

上述分类器识别方法的缺点为：这种分类方法对噪声点比较敏感，而且分割结果存在的误分类情况严重。

技术实现要素：

本发明的实施例提供了一种视网膜血管图像的分割方法，以实现有效地从视网膜血管图像中分割出粗血管段和细血管段。

为了实现上述目的，本发明采取了如下技术方案。

一种视网膜血管图像的分割方法，包括：

对视网膜血管图像进行双尺度匹配滤波处理，得到细尺度匹配滤波响应图像和粗尺度匹配滤波响应图像；

从所述细尺度匹配滤波响应图像中分割出线支持区域，使用局部自适应阈值方法对每一个线支持区域进行二值化处理，分割出细血管段；

应用固定比例阈值算法对所述粗尺度匹配滤波图像进行分割，得到粗血管段。

进一步地，所述的对视网膜血管图像进行双尺度匹配滤波处理，得到细尺度匹配滤波响应图像和粗尺度匹配滤波响应图像，包括：

提取彩色视网膜血管图像的RGB三个通道中的绿色通道；

用高斯函数来模拟视网膜血管的横切面灰度曲线，得到如下匹配滤波器：

式中，K(x,y)被称为核函数，σ是高斯函数沿x轴坐标中心的偏离度，L是高斯函数沿y轴被截断的闪电通道长度，式中x,y需满足|x|≤3σ,|y|≤L/2

以15°为间隔，选取角度区间[0°,180°]中的12个方向，创建12个匹配滤波器；

将所述彩色视网膜血管图像中的绿色通道分别与所述12个匹配滤波器做卷积计算，得到匹配滤波响应图像，将所述匹配滤波响应图像归一化并量化为256级的灰度图，当所述偏离度σ小于设定的阈值时，将得到的灰度图作为细尺度匹配滤波响应图像；当所述偏离度σ不小于设定的阈值时，将得到的灰度图作为粗尺度匹配滤波响应图像。

进一步地，所述的从所述细尺度匹配滤波响应图像中分割出线支持区域，包括：

计算细尺度匹配滤波图像中的每个像素点的梯度幅值和梯度方向，将所有像素点按照其梯度幅值大小进行排序，选取具有最高梯度幅值的像素点作为种子点，将梯度幅值小于设定的梯度阈值的像素点排除线支持区域的构建过程；

基于所述种子点利用区域生长算法生成若干个线支持区域，每个线支持区域包括一个种子点，并且为一个与种子点具有相似梯度方向的像素集合，每个像素点包括两个状态：使用过和未使用。

进一步地，所述的基于所述种子点利用区域生长算法生成若干个线支持区域，包括：

从像素点的排序列表中选择一个未使用的像素作为种子点，将所述种子点的梯度方向作为要生成的所述种子点所在的线支持区域的初始角度θ_region，将所述种子点的邻域中未使用的且其梯度方向跟区域角度θ_region之间的误差在τ之间的像素点添加到所述线支持区域中，根据更新后的像素点更新计算所述线支持区域的角度；

重复执行上述处理过程，直到所述种子点的邻域中没有符合条件的像素点添加到所述线支持区域中，对所述线支持区域对应的最小外接矩形进行扩展。

进一步地，所述的使用局部自适应阈值方法对每一个线支持区域进行二值化处理，分割出细血管段，包括：

将细尺度匹配滤波图像分割出多个线支持区域后，使用局部自适应阈值方法应用Otsu算法对每一个线支持区域进行二值化处理，分割出前景和背景，分割出单个的细血管段；

所述Otsu方法搜索最优的阈值使得前景与背景之间的方差最大，设t为前景和背景的分割阈值，则计算前景像素的概率w_0t和平均灰度u_0t，背景像素的概率w_1t和平均灰度为u_1t，前景和背景之间的方差表示为：

g_t＝w_0t·(u_0t-u_t)²+w_1t·(u_1t-u_t)²＝w_0t·w_1t·(u_0t-u_1t)²,t∈[0,255]

其中u_t表示图像总平均灰度，t的取值范围为0-255，当方差g_t最大时，前景和背景差异最大，则对应的灰度t是最佳阈值。

进一步地，所述的应用固定比例阈值算法对所述粗尺度匹配滤波图像进行分割，得到粗血管段，包括：

应用固定比例阈值算法对所述粗尺度匹配滤波图像进行分割，得到粗血管图像，所述固定比例阈值算法的阈值由以下公式计算：

其中r是输入参数，表示预期的血管比例，Num是频次计算函数，Total表示像素总数；

在应用所述固定比例阈值算法时，先对所述粗尺度匹配滤波图像中的像素进行降序排序，搜索最优的阈值Tr，根据所述阈值Tr对所述粗尺度匹配滤波图像进行二值化处理，分割出前景和背景，分割出单个的粗血管段。

进一步地，所述的方法还包括：

应用逻辑或操作对所述粗血管段的分割结果和所述细血管段的分割结果进行融合，得到完整的粗血管段和细血管段的分割结果。

由上述本发明的实施例提供的技术方案可以看出，本发明实施例的方法通过ALT方法可以有效地从视网膜血管图像中分割出细血管，通过FRT方法可以有效地从视网膜血管图像中分割出完整的粗血管，融合ALT方法和FRT方法可以得到完整的视网膜血管分割结果，分割结果准确率高。

本发明附加的方面和优点将在下面的描述中部分给出，这些将从下面的描述中变得明显，或通过本发明的实践了解到。

附图说明

为了更清楚地说明本发明实施例的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。

图1为本发明实施例提供的一种基于血管支持区域的双尺度非线性阈值的视网膜血管图像的分割方法的实现原理示意图；

图2为本发明实施例提供的一种基于血管支持区域的双尺度非线性阈值的视网膜血管图像的分割方法的处理流程图；

图3为本发明实施例提供的一种矩形扩展图，其中图(a)表示区域增长算法得到的矩形，图(b)表示横坐标和纵坐标两个方向上矩形的扩展，图(c)表示扩展之后的矩形；

图4为本发明实施例提供的一种线支持区域灰度直方图；

图5为本发明实施例提供的一种粗细血管融合图，图(a)表示局部自适应阈值的检测结果实例，图(b)表示固定比例阈值算法的分割结果，图(c)表示最终的融合结果。

具体实施方式

下面详细描述本发明的实施方式，所述实施方式的示例在附图中示出，其中自始至终相同或类似的标号表示相同或类似的元件或具有相同或类似功能的元件。下面通过参考附图描述的实施方式是示例性的，仅用于解释本发明，而不能解释为对本发明的限制。

本技术领域技术人员可以理解，除非特意声明，这里使用的单数形式“一”、“一个”、“所述”和“该”也可包括复数形式。应该进一步理解的是，本发明的说明书中使用的措辞“包括”是指存在所述特征、整数、步骤、操作、元件和/或组件，但是并不排除存在或添加一个或多个其他特征、整数、步骤、操作、元件、组件和/或它们的组。应该理解，当我们称元件被“连接”或“耦接”到另一元件时，它可以直接连接或耦接到其他元件，或者也可以存在中间元件。此外，这里使用的“连接”或“耦接”可以包括无线连接或耦接。这里使用的措辞“和/或”包括一个或更多个相关联的列出项的任一单元和全部组合。

本技术领域技术人员可以理解，除非另外定义，这里使用的所有术语(包括技术术语和科学术语)具有与本发明所属领域中的普通技术人员的一般理解相同的意义。还应该理解的是，诸如通用字典中定义的那些术语应该被理解为具有与现有技术的上下文中的意义一致的意义，并且除非像这里一样定义，不会用理想化或过于正式的含义来解释。

为便于对本发明实施例的理解，下面将结合附图以几个具体实施例为例做进一步的解释说明，且各个实施例并不构成对本发明实施例的限定。

本发明实施例提出了一种基于血管支持区域的双尺度非线性阈值的视网膜血管图像的分割方法，该方法首先运用两个不同尺度的高斯滤波器对眼底视网膜血管图像进行预处理操作，得到两个不同尺度的匹配滤波响应图像。然后分别对这两个不同尺度的匹配滤波响应图像进行视网膜血管提取。最后对提取到的血管进行合并，得到最终的结果。具体步骤如下：

本发明实施例提出了一种基于血管支持区域的双尺度非线性阈值的视网膜血管图像的分割方法的实现原理示意图如图1所示，具体处理流程如图2所示，包括如下的处理步骤：

步骤S210、对视网膜血管图像进行双尺度匹配滤波(double-scale filtering,DSF)处理，得到细尺度匹配滤波响应(fine matched filtering response,FMFR)图像和粗匹配滤波响应(coarse matched filtering response,CMFR)图像。

不同尺度的高斯匹配滤波器对于细血管和粗血管的增强效果是不一样的。更具体的，小尺度滤波器有利于突出细血管，提取粗血管的骨架；而粗尺度滤波器有利于增强粗血管，模糊化末梢的细血管。因此，我们设计了细尺度和粗尺度两类高斯核函数，将它们分别作用于视网膜血管图像的绿色通道，得到细尺度匹配滤波响应图像和粗尺度匹配滤波响应图像。

在彩色视网膜血管图像的RGB(Red、Green、Blue)三个通道中，绿色通道分量图像的血管与背景对比度最高，更加有利于血管分割，因此本发明首先提取彩色视网膜血管图像的绿色通道。

在视网膜血管图像中，血管中心像素点亮度较小，两边的像素点亮度较大，视网膜血管的横截面灰度轮廓可以用高斯型曲线近似。因此高斯匹配滤波方法常用来提升图像对比度。假定视网膜血管为分段等宽的直线段，其长度为L，宽度为2σ，我们用高斯函数来模拟视网膜血管的横切面灰度曲线，从而得到如下匹配滤波器：

式中，K(x,y)被称为核函数，σ是高斯函数沿x轴坐标中心的偏离度，L是高斯函数沿y轴被截断的闪电通道长度，为了使得匹配更精准，式中x,y需满足|x|≤3σ,|y|≤L/2。根据实验结果，我们设定L＝7。

因为血管方向是任意的，我们以15°为间隔考虑角度区间[0°,180°]中的12个方向，创建12个匹配滤波器。

视网膜血管图像分别与这12个高斯核做卷积，匹配滤波响应图像中的每一个像素的值等于最大的卷积值。为了便于后续处理，匹配滤波响应图像被归一化和量化为256级的灰度图。

我们可以观察到视网膜血管图像中既包含视盘附件的粗血管，也包含末梢的细血管。在应用匹配滤波时，如果选择比较小数值范围的σ，这个比较小范围是上文提到的血偏离度管的宽度或者说高斯函数沿x轴坐标中心的偏离度，大概为1.3～1.6个像素，则滤波结果图像中的细血管更加容易得到加强，粗血管被腐蚀；相反，如果选择比较大数值范围的σ，这个比较大数值范围是2.0～2.4个像素，则粗血管得到加强，细血管被模糊化。

因此，本发明实施例提出了双尺度匹配滤波方法。细尺度匹配滤波器选择较小的σ，增强细血管，同时抑制噪声和平滑背景区域。细尺度匹配滤波器产生的响应结果经过量化，得到FMFR图像，该图像将是自适应阈值分割算法的输入。相反，粗尺度匹配滤波器选择较大的σ，增强粗血管部分，得到CMFR图像。因为粗血管的边缘部分在细尺度匹配滤波图像中容易被腐蚀，但是粗血管比较容易从粗尺度匹配滤波图像中完整分割出来，所以，CMFR图像将用于分割结果融合。

步骤S220、从细尺度匹配滤波响应图像中分割出线支持区域(vesselsupport region，VSR)，使用局部自适应阈值(adaptive local thresholding,ALT)方法应用Otsu算法对每一个VSR进行二值化处理，分割出单个的细血管段。

经过匹配滤波后，FMFR图像的对比度得到了增强，尤其是污点和病变区域得到了抑制。但是FMFR图像中血管的灰度分布还是比较分散，部分血管的灰度值与背景灰度值存在较大重叠。在理论上，我们找不到一个全局阈值线性分割血管和背景。VSR是指包含一个血管段的矩形区域，可以通过算法自动检测。在一个局部VSR中，其直方图具有明显的双模态性质，我们可以应用自动阈值算法(比如Otsu)分割血管和背景。该过程首先自动检测FMFR中的所有VSR区域，然后应用Otsu算法分割每一个VSR区域，同时把所有非VSR区域的像素设置为背景，最后得到细血管分割图(fine vessel segmentat ion,FVS)。

S2-1:线支持区域检测

由于视网膜细血管对小尺度参数σ比较敏感，所以本发明首先计算在尺度σ＝1.3条件下的匹配滤波图像，然后在此基础上进行以下处理。

S2-1-1：线支持区域生成

首先计算细尺度匹配滤波图像中的每个像素点的梯度幅值和梯度方向，然后将所有像素点按照其梯度幅值大小进行排序。较强的边缘点或区域一般都具有比较高的梯度幅值，通常在视网膜血管边缘的像素具有最高的梯度幅值，因此首先选取具有最高梯度幅值的像素点作为种子点。在计算过程中，梯度幅值小于q(本发明选用0.4)的像素点将被拒绝参与线支持区域的构建过程。

本发明利用区域生长算法生成若干个线支持区域，每个线支持区域也即一个与种子点具有相似梯度方向的像素集合。每个像素点包括两个状态，即使用过和未使用。初始状态将所有像素点全部置为未使用。

区域生长算法首先从排序列表中选择一个未使用的像素作为种子点，该像素的邻域中未使用的且其梯度方向跟区域角度θ_region之间的误差在τ之间的像素将被加入到该区域中。文中的试验τ的范围大概在18°到24°之间，一般默认取22。

区域的初始角度θ_region就是种子点的梯度方向，每次添加一个新的像素到该区域，区域的角度就被更新。区域的角度就被更新为：

θ_j表示像素点j梯度的垂直方向。

i(x,y)表示像素(x,y)点处的灰度值，g_x(x,y)、g_y(x,y)分别表示像素(x,y)在x和y方向的梯度值。

如此依次进行，直到没有任何像素可以添加到矩形当中。

前面得到的线支持区域，用一个最小外接矩形来表示。从而可以获取矩形的一些基本信息，图3为本发明实施例提供的一种矩形扩展图，包括矩形中心点的坐标以及矩形的长度和宽度以及主方向，其中图(a)表示区域增长算法得到的矩形，图(b)表示横坐标和纵坐标两个方向上矩形的扩展，图(c)表示扩展之后的矩形。

矩形的扩展，这里分为两步，第一步：先对矩形横向和纵向都进行等值扩展，扩展幅度选取为矩形宽度width的一半)。这样每次扩展完一个矩形后，将矩形里的像素点的状态都设置为Used，下次这些被设置为Used的点就不会被选作种子点。第二步：再对之前的扩展的矩形基础上，对横向和纵向都进行等值扩展，扩展幅度和之前的定值大小一样。在扩展后的矩形基础上进行阈值处理，这样就会解决两个矩形之间没有交叉的问题，因为在第二次扩展之后增加的那些像素可以被选作种子点。

算法1.区域生长算法

S2-1-2:矩形生长

区域生长算法得到的线支持区域能很好地覆盖大部分闪电通道区域，但是它存在两方面的不足：(1)闪电通道两侧的梯度值比较大，因而区域增长的种子点一般为闪电通道两个边缘的像素点，导致在闪电通道同一横截面处会出现两个矩形；(2)下一个区域增长的种子点不一定是从已生成的矩形区域中选择的，因而新生成的区域跟上一个区域不一定有重叠区域。也即沿着闪电通道方向矩形可能不连续，从而不能将所有的前景区域全部包含进去。

因此，本发明提出了矩形扩展算法，即在原来的矩形基础上进行双向扩展，扩展幅度选取为矩形宽度的一半。其中纵向扩展(沿着血管发展的方向)可以解决矩形区域不连续的问题，横向扩展(与血管发展方向垂直)可以达到合并两个矩形的效果。

图3给出了矩形扩展的示意图，其中O点为矩形中心点，theta为矩形的主方向，Length和Width分别为矩形的长和宽，P、Q分别为矩形宽上的中心点，可根据x1和x2,以及y1和y2的大小分为9种情况。

S2-2:线支持区域阈值分割

应用上述VSR检测算法，一幅MFR图像可以分割出多个VSR。

步骤S230、然后ALT方法应用Ot su算法对每一个VSR进行二值化处理，分割出前景和背景，分割出单个的细血管段。直观的，每一个VSR区域都是对比度非常明显的图像块，血管段与背景在灰度空间具有明显区别。从统计角度分析，VSR区域的强度值分布具有双模态性，即背景像素和血管像素集中分布在各自的强度区间。图4展示了随机抽取的两个VSR的灰度直方图。我们的更多实验结果都表明VSR的强度分布具有双模态性。

Otsu是一种经典的自动阈值方法，它对具有双模态分布的图像具有较好的分割效果。Otsu方法的原理是搜索最优的阈值使得前景与背景之间的方差最大。假设t为前景和背景的分割阈值，则可以计算前景像素的概率w_0t和平均灰度u_0t，背景像素的概率w_1t和平均灰度为u_1t。前景和背景之间的方差可表示为：

g_t＝w_0t·(u_0t-u_t)²+w_1t·(u_1t-u_t)²＝w_0t·w_1t·(u_0t-u_1t)²,t∈[0,255]其中u_t表示图像总平均灰度，t的取值范围为0-255。当方差g_t最大时，前景和背景差异最大，则对应的灰度t是最佳阈值。

步骤S240、应用固定比例阈值算法分割CMFR图像，得到粗血管分割图，融合粗细血管分割方法，得到完整的细血管和粗血管的分割结果。

融合包括两个主要步骤：首先，应用固定比例阈值算法分割CMFR图像，得到粗血管分割图(coarse vessel segmentation,CVS)。然后，FVS和CVS通过逻辑或操作进行融合，使得融合结果既包含细血管，也包含完整的粗血管。

ALT可以检测得到细血管图像，但是ALT分割的粗血管往往只包含其骨架，而遗漏了其外围部分的像素。图5(a)展示了ALT的检测结果实例。为了提高ALT的检测性能，本发明提出细血管和粗血管的融合方法。该方法首先应用固定比例阈值算法(Fixed-ratio thresholding，FRT)对粗尺度匹配滤波图像进行分割，得到粗血管图像，然后融合细血管图像和粗血管图像，得到最终的粗细血管融合结果。

固定比例阈值算法是一种简单的基于先验的二值化方法。直观的，视网膜图像具有明显的结构性，即它由线状的血管和平坦的背景组成，而且血管部分的比例往往比较低。从统计角度分析，DRIVE中血管像素比例的平均值和方差分别为8.43％和1.38％，STARES中血管像素比例的平均值和方差分别为7.6％和3.15％。2个数据集DRIVE和STARES都只用平均值作为评价指标，去掉方差，DRIVE的平均值12.7％，STARES的平均值10.4％。

需要提出的是，STARES中有一些病变的视网膜图像，所以其血管比例的方差比较大。在粗尺度匹配滤波图像中，粗血管像素的响应值比细血管和背景像素的响应值大。因此，FRT方法的阈值可以由以下公式计算：

其中r是输入参数，表示预期的血管比例。Num是频次计算函数，Total表示像素总数。在实现该算法时，先应用桶排序算法对粗尺度匹配滤波图像中像素进行降序排序，然后搜索最优的阈值Tr，最后根据Tr对粗尺度匹配滤波图像进行二值化。图5(b)展示FRT的分割结果，可以观察到该结果较完整的分割出粗血管，尽管它丢失了细血管的细节。

综上，ALT善于分割细血管，FRT则可以分割出完整的粗血管。因此，融合两个方法的结果可以预期比较完善的分割性能。简单的，本发明应用逻辑或操作对于ALT和FRT的结果进行融合。即两幅图像中对应位置像素点的灰度值有一个为255，那么结果图像中对应位置像素点的灰度值即为255，只有当两幅图像中对应位置像素点的灰度值均为0，结果图像中对应位置像素点的灰度值才为0。

另外，为了消除背景噪声和部分病变组织的干扰，去除面积小于10个像素点的区域。图5(c)展示了最终的融合结果。

综上所述，本发明实施例的方法通过ALT方法可以有效地从视网膜血管图像中分割出细血管，通过FRT方法可以有效地从视网膜血管图像中分割出完整的粗血管，融合ALT方法和FRT方法可以得到完整的视网膜血管分割结果，分割结果准确率高。

目前国内外大部分方法都只针对正常的、成像较好的视网膜图像进行血管提取，而对于低对比度或者发生病变的视网膜图像，由于血管和背景区域像素灰度值大小较为接近，已有的传统方法大部分无法将血管与背景正确地分割出来。而本发明利用大小尺度的高斯匹配滤波将血管分为粗细血管，分别进行增强处理，效果明显。对于细血管，利用了像素点的梯度大小和方向，基于区域增长和矩形扩展的方法，具有较好的局部自适应性，可以快速地判定出前景区域。而对于粗血管，使用全局固定阈值，可以较好地保留主干部分。

本领域普通技术人员可以理解：附图只是一个实施例的示意图，附图中的模块或流程并不一定是实施本发明所必须的。

通过以上的实施方式的描述可知，本领域的技术人员可以清楚地了解到本发明可借助软件加必需的通用硬件平台的方式来实现。基于这样的理解，本发明的技术方案本质上或者说对现有技术做出贡献的部分可以以软件产品的形式体现出来，该计算机软件产品可以存储在存储介质中，如ROM/RAM、磁碟、光盘等，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)执行本发明各个实施例或者实施例的某些部分所述的方法。

本说明书中的各个实施例均采用递进的方式描述，各个实施例之间相同相似的部分互相参见即可，每个实施例重点说明的都是与其他实施例的不同之处。尤其，对于装置或系统实施例而言，由于其基本相似于方法实施例，所以描述得比较简单，相关之处参见方法实施例的部分说明即可。以上所描述的装置及系统实施例仅仅是示意性的，其中所述作为分离部件说明的单元可以是或者也可以不是物理上分开的，作为单元显示的部件可以是或者也可以不是物理单元，即可以位于一个地方，或者也可以分布到多个网络单元上。可以根据实际的需要选择其中的部分或者全部模块来实现本实施例方案的目的。本领域普通技术人员在不付出创造性劳动的情况下，即可以理解并实施。

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