基于动态尺度分配的视网膜血管提取方法及系统与流程

本发明涉及基于动态尺度分配的视网膜血管提取方法及系统。

背景技术：

迄今为止，常用的视网膜血管自动提取算法有：

1、基于视网膜血管追踪方式的算法。这类方法可以较为完整的提取出视网膜血管的网络，但是算法复杂度较高，运算量较大。此外，对于一些对比度较低的视网膜血管图像，这类算法的提取准确度不够。其中比较典型的视网膜血管追踪算法是由Tolias在1998年提出的基于模糊C均值聚类算法，是在血管的起始处(视盘)挑选出合适的种子点，并由此对整个视网膜血管网络进行追踪。建立起视网膜血管横截面的一维模型，并通过建立起种子点与血管一维模型之间的模糊相似关系，对所有种子点进行分类，判断其是否属于血管，从而完成最终的血管网络分割提取。此类算法存在的最大弊病是种子点的选取会直接影响提取结果，而且在视盘中选取一个最合适的像素点作为种子点是一个比较繁重的工作。此外，此类算法对于血管的分支点的处理效果欠佳，并且会难以避免地丢失细小血管结构，使得分割结果不够准确。

2、基于分类器的提取算法。这类方法的主要思想是利用视网膜血管所提供的一些先验信息，构造出一个合适的分类模型，即通常意义上的分类器，并利用构造好的分类器将视网膜图像中的像素点进行映射分类。Boyce(1999)介绍了一种监督方式的视网膜血管提取算法。这种基于分类器的识别方式对噪声极为敏感，所以最终的分类效果不是很好。

3、基于匹配滤波器的提取算法，是最为经典也是使用最为广泛的一种方式。根据血管的灰度分布特征，大多数匹配滤波算法均选择高斯滤波器与视网膜图像进行匹配滤波，生成的响应具有较高的输出信噪比。其中，最早使用二维高斯匹配滤波器对视网膜血管结构进行提取的方法是Chaudhuri et al.在1989年提出的。此算法基于一个事实：视网膜血管横截面的灰度分布是服从高斯分布的。因此，若将建立好的高斯滤波器模板与视网膜血管图像进行不同方向的匹配滤波，当滤波器的尺度与血管的宽度在一定范围内相吻合时，相应宽度的血管就得以增强，产生很大的卷积响应。但是文献中仅采用了一种尺度的匹配滤波器，无法使所有宽度的血管均得以增强，并且不论是血管结构还是非血管结构均对高斯匹配滤波器有很强的响应，这样会导致错误分割像素点的出现。针对经典算法的弊端，多种多尺度匹配滤波方案相继被提出，这一改进旨在利用不同尺度的滤波器对不同宽度的血管进行增强，并在后续步骤中将其准确提取。Bob Zhang在2010年提出的基于高斯一阶导数滤波器的视网膜血管提取算法(MF-FDOG)，是在Chaudhuri et al.所提出算法的基础上，结合高斯一阶导数型匹配滤波器与图像产生卷积响应的特点，使用了多尺度的匹配滤波器对所有视网膜图像像素点进行遍历分类。这种方法虽然在一定程度上克服了台阶边缘对分割结果造成的影响，但是分割结果中仍有噪声模式存在，一些被噪声淹没的细小血管也未被提取出来。此外，此方法对于血管分支点与连通性的处理也有待提高。Qin Li在2012年提出了基于多尺度滤波响应乘积(MPMF)的方法。此方法利用各尺度响应的乘积将各个宽度的血管信息在尺度域进行混合，不仅对血管进行了有效的增强且抑制了噪声，同时对分割的血管有较好的宽度估计。但本方法仍然存在一定的局限性：其使用了三个不同尺度的滤波器对视网膜血管进行提取，对于变化较大的血管宽度而言，仍存在某些不能被很好增强的血管宽度，因此难以准确地分割出所有血管；另外对于一些存在病变区域视网膜图像，还会有错误分割像素点的出现。

欲以现有方法得到更高的视网膜血管分割精度，目前难以解决的技术问题是：

1、如何利用多尺度方案使各宽度的血管尽量完整地被提取出来；

2、如何减少错误分割点以及抑制噪声；

3、如何准确地估计血管宽度。

技术实现要素：

本发明的目的就是为了解决上述问题，提供基于动态尺度分配的视网膜血管提取方法及系统，通过对视网膜图像的预处理、图像分块、血管分类、动态尺度分配、多尺度匹配滤波、阈值分割和后处理，实现对视网膜图像的血管提取，并在剔除复杂非血管结构的同时，避免了对血管宽度的过分估计，实现了更简单、更准确的视网膜血管提取。

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

基于动态尺度分配的视网膜血管提取方法，包括如下步骤：

步骤(1)：视网膜图像预处理：对彩色视网膜图像的绿色通道分量进行对比度增强；

步骤(2)：图像分块：对预处理后的视网膜图像分割成设定个数的子图像；

步骤(3)：血管分类：将每个子图像中血管分为大中小三类；

步骤(4)：动态尺度分配：动态选择不同尺度的滤波器对不同宽度血管进行增强；

步骤(5)：多尺度匹配滤波：采用多尺度方案并利用高斯匹配滤波模板和高斯一阶导数匹配滤波模板对视网膜血管图像进行滤波处理；

步骤(6)：阈值处理：提取出血管结构并剔除非血管结构，将所有子图像的提取结果进行重新拼接，得到视网膜血管网络二值图像；

步骤(7)：后处理：对阈值处理后的视网膜血管网络二值图像进一步消除噪声，对血管边缘进行平滑处理，消除图像中残留的视网膜边界，对不连续的细小血管进行断点连接，从而保留图像细节，得到分割精度高的视网膜血管网络图像。

进一步的，所述步骤(1)的步骤为：

步骤(1-1)：提取彩色视网膜图像的绿色通道分量：彩色视网膜图像含有红色、绿色和蓝色三个通道，仅选择对比度高、噪声低的绿色通道作为初始处理对象；

步骤(1-2)：多尺度顶帽变换：利用形状不变，尺寸等差增大的圆形结构元素，对初始处理对象进行顶帽变换处理，增强初始处理对象的对比度；

步骤(1-3)：基于高斯曲线拟合的直方图线性拉伸：对经多尺度顶帽变换增强所得图像的灰度直方图进行基于高斯曲线拟合的直方图线性拉伸，得到预处理后的视网膜图像。

进一步的，所述步骤(2)的步骤为：

对预处理后的视网膜图像按照尺寸进行均等分块，分割成设定个数的子图像。

进一步的，所述步骤(3)的步骤为：

根据血管宽度将血管分为大中小三类，依照不同宽度血管的灰度分布特征，对子图像包含的血管种类进行判断。

进一步的，步骤(3)对步骤(2)的子图像自上而下、自左向右进行编号；先根据子图像位置的不同将子图像分为三类：位于原图四个角的子图像、位于原图中央的子图像和位于原图四条边的子图像；然后针对每类子图像，根据不同宽度血管的灰度分布特征，将每个子图像的血管进行分类。

进一步的，针对每类子图像，根据不同宽度血管的灰度分布特征，将每个子图像的血管进行分类，步骤为：

(3-1)对于位于原图四个角的子图像，子图像中均仅含有中血管；

(3-2)对于位于原图中央的子图像，

判断是否存在大血管，对位于原图中央的子图像灰度直方图的第一设定灰度分布范围内像素点出现的频率进行统计，当统计值大于第一设定阈值时，则代表该子图像中含有大血管，反之，不含大血管；

判断是否存在中小血管，对位于原图中央的子图像灰度直方图的第二设定灰度分布范围内的灰度分布进行高斯曲线拟合，得到高斯函数峰值和方差：

若高斯函数峰值大于第二设定阈值且方差小于第三设定阈值，则子图像第二设定灰度分布范围内有小血管但无中血管；

若高斯函数峰值小于第二设定阈值且方差小于第三设定阈值，则子图像第二设定灰度分布范围内既无小血管又无中血管；

若高斯函数峰值大于第二设定阈值且方差大于第三设定阈值，则子图像第二设定灰度分布范围内既有小血管又有中血管；

若高斯函数峰值小于第二设定阈值且方差大于第三设定阈值，则子图像第二设定灰度分布范围内无小血管但有中血管。

(3-3)对于位于原图四条边的子图像，

判断是否存在大血管，对位于原图边缘的子图像灰度直方图的第一设定灰度分布范围内像素点出现的频率进行统计，当统计值大于第四设定阈值时，则代表该子图像中含有大血管，反之，不含有大血管；

判断是否存在中小血管的标准，与(3-2)类似，但阈值设定不同，对位于原图边缘的子图像灰度直方图的第二设定灰度分布范围内的灰度分布进行高斯曲线拟合，得到高斯函数峰值和方差：

若高斯函数峰值大于第五设定阈值且方差小于第六设定阈值，则子图像第二设定灰度分布范围内有小血管但无中血管；

若高斯函数峰值小于第五设定阈值且方差小于第六设定阈值，则子图像第二设定灰度分布范围内既无小血管又无中血管；

若高斯函数峰值大于第五设定阈值且方差大于第六设定阈值，则子图像第二设定灰度分布范围内既有小血管又有中血管；

若高斯函数峰值小于第五设定阈值且方差大于第六设定阈值，则子图像第二设定灰度分布范围内无小血管但有中血管。

所述步骤(4)的步骤为：

利用滤波器尺度与血管宽度之间的关系，动态选择不同尺度的滤波器对不同宽度血管进行增强。

进一步的，所述步骤(4)的步骤为：

步骤(4-1)：提取大血管所采用的滤波器尺度参数为2；

步骤(4-2)：提取中血管所采用的滤波器尺度参数为1和0.7，提取小血管所采用的滤波器尺度参数为0.5和0.2。

所述步骤(5)的步骤为：

采用高斯匹配滤波模板和高斯一阶导数匹配滤波模板同时对增强后的视网膜图像进行匹配滤波，分别得到高斯匹配滤波器与视网膜图像的卷积响应及高斯一阶导数匹配滤波器与视网膜图像的卷积响应；并进一步计算得到高斯一阶导数匹配滤波器与视网膜图像卷积响应的局部均值响应。

所述步骤(5)的高斯匹配滤波模板：

高斯型匹配滤波模板上坐标为(x,y)的一点，其权重表示如下：

其中，σ₁表示高斯匹配滤波器的尺度，L代表与滤波器模板瞬时方向相平行的血管长度，m代表高斯滤波器模板中系数的均值，

式中Q表示构建的高斯滤波器模板中所包含的点数，N表示滤波器模板所在的邻域，表示高斯滤波器模板中的一个离散点。

在高斯型匹配滤波模板中减去m项，实现对视网膜图像背景的平滑处理。并且，由于视网膜血管在平面中的各个方向上均有分布，故而在匹配滤波过程中，滤波器模板需要不断的旋转以检测不同方向上的血管。

所述步骤(5)的高斯一阶导数匹配滤波模板：

高斯一阶导数匹配滤波模板上坐标为(x,y)的一点，其权重表示如下：

其中，σ₂表示高斯一阶导数匹配滤波器的尺度，L代表与滤波器模板瞬时方向相平行的血管长度。

所述步骤(6)的步骤为：

利用高斯一阶导数匹配滤波器与视网膜图像卷积响应的局部均值响应对每个像素点的分割阈值进行调整，得到整个图像的阈值矩阵；将整个图像的阈值矩阵与高斯匹配滤波器和视网膜图像的卷积响应进行逐点比较，提取出血管结构并剔除非血管结构，将所有子图像的提取结果进行重新拼接，得到视网膜血管网络二值图像；

进一步的，所述步骤(6)的步骤为：

定义A为高斯匹配滤波器与视网膜图像的卷积响应，为高斯一阶导数匹配滤波器与视网膜图像卷积响应的局部均值响应；

将所有子图像利用其灰度分布特征，按照血管与背景的灰度对比度的高低，以及是否含具有台阶边缘的非血管结构进行分类，对各子图像进行动态的阈值处理；

通过将整个图像的阈值矩阵和高斯匹配滤波器与视网膜图像的卷积响应逐点比较，每个像素点均与对应点的阈值进行比较，大于阈值的像素点被分类为血管，否则被分类为非血管；

最后，将各子图像的不同尺度分割结果利用逻辑“或”结合起来，将所有子图像的分割结果进行重新拼接，得到阈值处理后的视网膜血管网络二值图像。

所述将所有子图像利用其灰度分布特征，按照血管与背景的灰度对比度的高低，以及是否含具有台阶边缘的非血管结构进行分类的步骤为：

步骤(6-1)：利用步骤(3)血管分类过程中得到的高斯曲线方差作为判断图像对比度高低的依据；

方差c小于设定值c_th，表示子图像血管与背景的灰度对比度高，血管提取容易；

方差c大于设定值c_th，表示子图像血管与背景的灰度对比度低，血管提取困难。

对低对比度子图像需进行如下处理：对第二设定灰度分布范围进行分段处理，以像素点出现频率最高的灰度值为界，分为高低两个灰度级范围。

步骤(6-2)：

对图像灰度直方图第三设定灰度分布范围内像素点出现的频率和s进行统计，并设立第七设定阈值，频率大于第七设定阈值时，代表子图像中含具有台阶边缘的非血管结构。

若判断出含具有台阶边缘的非血管结构，那么接下来要对其进行定位。进一步对局部均值响应的各点进行遍历判断，当中坐标为(x,y)的某一任意点响应幅度大于b_th时，认为(x,y)位置附近将会有具有台阶边缘的非血管结构像素点出现；小于b_th时则认为(x,y)位置附近将会有血管像素点出现。提取步骤(3)中所判定的三种不同规格血管时，所设置的b_th值也是不同的：

b_th＝1小血管

b_th＝1.1中血管

b_th＝1.2大血管

步骤(6-3)：对各子图像的分割阈值由多阈值公式确定：

其中，T_k,j为子图像的阈值矩阵，h_i,j为常系数，为高斯滤波响应A的局部均值响应，k为图像血管与背景的灰度对比度对阈值个数的影响，j为有无台阶边缘对阈值个数的影响，K(c,c_th)与J(s,s_th)的取值具体如下：

其中，c是高斯拟合曲线的方差，当c小于阈值c_th时,表明对比度高，在整个灰度级范围内选择阈值，此时K的取值为1；反之，在对比度低的情况下，分为高低两个灰度级范围选择阈值，此时K的取值为2。s是第三设定灰度分布范围内像素点出现的频率和，如果s小于s_th,则不存在台阶边缘，不需对阈值针对血管结构和非血管结构进行区分，此时J的取值为1；当s大于等于s_th时,存在台阶边缘，血管结构和非血管结构要分别设置不同的阈值，此时J的取值为2。

综上所述，针对不同子图像，有以下四类阈值处理手段：

第一类：对于对比度低且存在台阶边缘的子图像，采用h₁₁、h₁₂、h₂₁、h₂₂四个阈值，其中h₁₁针对高灰度级范围内的血管结构，h₁₂针对高灰度级范围内的非血管结构，h₂₁针对低灰度级范围内的血管结构，h₂₂针对低灰度级范围内的非血管结构；

第二类：对于对比度低且不存在台阶边缘的子图像，采用h₁₁、h₂₁两个阈值，其中h₁₁针对高灰度级范围内的血管结构，h₂₁针对低灰度级范围内的血管结构；

第三类：对于对比度高且存在台阶边缘的子图像，采用h₁₁、h₁₂两个阈值，其中h₁₁针对第二设定灰度分布范围内的血管结构，h₁₂针对第二设定灰度分布范围内的非血管结构；

第四类：对比度高且不存在台阶边缘，采用h₁₁一个阈值，h₁₁针对整个灰度级范围内的血管结构。

整个图像的阈值矩阵为最终可以通过公式(7)得到分割后的血管网络：

其中，vess表示各像素点的最终取值。即将整幅图像的阈值矩阵与图像的高斯滤波响应A逐点比较，每个像素点均利用对应点的阈值判断一次，大于对应点阈值的像素点被分类为血管，否则被分类为非血管。

步骤(6-4)：将各子图像的不同尺度分割结果利用逻辑“或”操作结合起来，再将所有子图像的分割结果进行重新拼接，得到阈值处理后的视网膜血管网络二值图像。

所述步骤(7)的步骤为：

首先，利用目标视网膜血管与噪点所在区域的几何特征差别，先删除视网膜血管网络二值图像中面积小于设定阈值P的连通域，再利用连通域外接矩形的长宽比，进一步的消除图像中残留的噪声；

其次，使用多尺度的高斯滤波器模版来对阈值分割后的血管边缘进行平滑处理；

接下来，利用DRIVE数据库中提供的掩膜图像，找到视网膜边界像素点的位置，进而消除图像中残留的视网膜边界；

最后，采用闭运算对视网膜血管网络二值图像中不连续的细小血管进行断点连接。

基于动态尺度分配的视网膜血管提取系统，包括：

视网膜图像预处理模块：对彩色视网膜图像的绿色通道分量进行对比度增强；

图像分块模块：对预处理后的视网膜图像按照尺寸进行均等分块，分割成设定个数的子图像；

血管分类模块：根据血管宽度将血管分为大中小三类。并依照不同宽度血管的灰度分布特征，对子图像包含的血管种类进行判断。

动态尺度分配模块：利用滤波器尺度与血管宽度之间的关系，动态选择不同尺度的滤波器对不同宽度血管进行增强；

多尺度匹配滤波模块：采用高斯匹配滤波模板和高斯一阶导数匹配滤波模板同时对增强后的视网膜图像进行匹配滤波，分别得到高斯匹配滤波器与视网膜图像的卷积响应及高斯一阶导数匹配滤波器与视网膜图像的卷积响应；并进一步计算得到高斯一阶导数匹配滤波器与视网膜图像卷积响应的局部均值响应；

阈值处理模块：利用高斯一阶导数匹配滤波器与视网膜图像卷积响应的局部均值响应对每个像素点的分割阈值进行调整，得到整个图像的阈值矩阵。将整个图像的阈值矩阵与高斯匹配滤波器和视网膜图像的卷积响应进行逐点比较，提取出血管结构并剔除非血管结构，将所有子图像的提取结果进行重新拼接，得到视网膜血管网络二值图像；

后处理模块：对阈值处理后的视网膜血管网络二值图像进一步消除噪声，对血管边缘进行平滑处理，消除图像中残留的视网膜边界，对不连续的细小血管进行断点连接，从而保留图像细节，得到分割精度高的视网膜血管网络图像。

本发明的有益效果：

1、如果只利用高斯滤波器与视网膜图像所得到的卷积响应进行阈值分割，某些非血管结构就会不可避免地被错误分割出来，致使最终的分割精度受到较大影响。因此，本发明使用两个匹配滤波模版同时作用，来减小非血管结构被错误分割的可能性。

2、匹配滤波部分是在对基于高斯一阶导数滤波器的视网膜血管提取算法(MF-FDOG)与基于多尺度滤波响应乘积的视网膜血管提取算法(MPMF)的基础上，将整个视网膜图像进行分割处理，并针对各子图像的灰度分布特征动态的为其分配不同的滤波器尺度。

3、在多尺度方面，为了避免对中小血管的种类产生错误的判断，本发明分别为提取中小血管分配了两个不同的尺度参数。本发明一共提供了五个尺度参数，分别用于对大、中、小血管进行增强，是目前多尺度方案中设置尺度数最多的算法。使用这种动态分配尺度的多尺度滤波算法，在抑制具有台阶边缘的非血管结构以及保证血管宽度的同时，使得更多被噪声淹没的细小血管得以增强，进而在后续处理中将其分割提取出来，提高了算法的性能。

4、在阈值处理方面，本发明针对各子图像对比度的不同以及是否包含台阶边缘这几种情况进行分类处理。这种根据子图像不同的灰度分布情况动态的进行阈值处理的方法，使得对比度较低区域中的一些细小血管也能被分割提取出来。

5、后处理方面，本发明在消除图像噪声以及错误分割点的基础上，较好的保留了图像的细节，使血管分割精度得到了进一步的提高。

6、本发明算法具有较低的复杂度且易于实现，对推进视网膜疾病的诊断具有重大意义。

附图说明

图1为本发明的流程图；

图2为本发明的系统架构图；

图3(a)为DRIVE数据库测试集中的第19号彩色视网膜血管图像；

图3(b)为提取出的绿色通道图像；

图3(c)为多尺度顶帽变换后的图像；

图3(d)为灰度拉伸后的图像。

图4(a)为DRIVE数据库测试集中的第07号彩色视网膜血管图像；

图4(b)和图4(c)为图像分割后的两个子图像；

图5为各子图像的编号情况；

图6(a)为位于原图四角的子图像之一；

图6(b)为位于原图四角的子图像之二；

图7(a₁)为中央子图像a；

图7(a₂)为中央子图像a的灰度分布直方图；

图7(b₁)为中央子图像b；

图7(b₂)为中央子图像b的灰度分布直方图；

图7(c₁)为中央子图像c；

图7(c₂)为中央子图像c的灰度分布直方图；

图7(d₁)为中央子图像d；

图7(d₂)为中央子图像d的灰度分布直方图；

图8为利用高斯拟合曲线的峰值a和方差c来对中小血管划分示意图；

图9(a)为DRIVE数据库测试集中的第19号彩色视网膜血管图像；

图9(b)为预处理后的图像；

图9(c)为本发明算法提取的血管网络；

图9(d)为ground truth血管网络图。

具体实施方式

下面结合附图与实施例对本发明作进一步说明。

如图1所示，基于动态尺度分配的视网膜血管提取方法，包括如下步骤：

步骤(1)：图像预处理

进一步的，所述步骤(1)的步骤为：

步骤(1-1)：提取彩色图像的绿色通道分量：原始彩色视网膜图像含有红色、绿色和蓝色三个通道，仅选择对比度高、噪声低的绿色通道作为初始处理对象；

步骤(1-2)：多尺度顶帽变换：利用形状不变，尺寸等差增大的圆形结构元素，对初始处理对象进行多尺度顶帽变换处理，增强初始处理对象的对比度；

步骤(1-3)：基于高斯曲线拟合的直方图线性拉伸：对经多尺度顶帽变换增强所得图像的灰度直方图进行基于高斯曲线拟合的直方图线性拉伸，得到预处理后的视网膜图像。

通过预处理步骤增强后的视网膜图像如图3(a)、图3(b)、图3(c)和图3(d)所示。

步骤(2)：图像分块：对预处理后的视网膜图像按照尺寸进行均等分块，分割成设定个数的子图像；

由于视网膜图像中血管和背景在局部范围内都比较均匀，所以将整个视网膜图像均等的分割成20个等尺寸的子图像。图4(a)、图4(b)和图4(c)为增强后的视网膜图像及其子图像。

步骤(3)：根据血管宽度将血管分为大中小三类。并依照不同宽度血管的灰度分布特征，对子图像包含的血管种类进行判断。

将图像等分为20个无重叠的子图，并对分块自上而下、自左向右进行编号，如图5所示。根据位置的不同将20个子图像分为三类；针对血管宽度的不同，将血管分为大中小三类，并在每一个子图中分别判断其中的血管种类。

第一类：位于原图四角的子图像，编号1、4、17、20。少部分属于这一类的子图像图6(a)所示，其中含有对比度较高的血管，但属于血管的像素所占比重较小；多数子图像均如图6(b)所示，其中所包含的血管对比度极低，甚至没有血管像素，所以很难通过分析其直方图的灰度分布特征对所含血管的种类进行判断。针对数据库中所有图像进行分析，且为了权衡上述两种情况，本算法统一为其分配中等规格的尺度。

第二类：处于原图像中央的子图像，编号为6、7、10、11、14、15。

此类子图像处于原图像中央，因而其中不包含视网膜之外的黑色背景；且图中的大血管较暗，所以属于大血管的像素点应处于直方图的低灰度级范围内。

1)判断是否存在大血管。对直方图第一设定灰度分布范围内像素点出现的频率进行统计，当统计值大于第一设定阈值时，则代表该子图像中含有大血管，反之，不含大血管；图7(a₁)和图7(b₁)分别代表了两种情况：图7(a₁)低灰度级范围内存在大血管，图7(b₁)低灰度级范围内不存在大血管。通过对其相应的灰度直方图图7(a₂)、图7(b₂)进行对比可以发现：含有大血管的子图像，在直方图第一设定灰度分布范围内像素点出现的频率明显高于不含有大血管的子图像。

2)判断是否存在中小血管。图7(c₁)和图7(d₁)分别代表了两种情况：图7(c₁)高灰度级范围内同时存在中、小血管，图7(d₁)高灰度级范围内仅存在小血管。通过对其相应的灰度直方图图7(c₂)、图7(d₂)进行对比可以发现：灰色背景以及中小血管的灰度分布范围基本集中在第二设定灰度分布范围之间，且其中小血管与灰色背景的灰度分布范围较为接近，甚至部分重叠。此时，对直方图高灰度级范围内的灰度分布进行高斯曲线拟合，得到的高斯函数峰值a和方差c可以作为区分中小血管的标准。若峰值大，则证明此时子图像的高灰度级范围内不仅存在灰色背景还存在小血管，使得在二者灰度分布重叠的范围内像素点出现的概率较大；若峰值小，则证明此时子图像的高灰度级范围内仅存在灰色背景；方差大时，证明在高灰度级范围内像素点在各灰度级上分布较为分散。因为中血管与小血管或灰色背景的灰度分布范围重叠较小，所以此时可以认为该子图像的高灰度级范围内存在中血管；方差小时，则认为该子图像的高灰度级范围内不存在中血管。图7(c₂)和图7(d₂)使上述区分标准得以验证。

根据以上方差与峰值大小不同的取值，对于子图像中是否含有中小血管的判断可以分为以下四类，如图8所示：

若高斯函数峰值大于第二设定阈值且方差小于第三设定阈值，则子图像第二设定灰度分布范围内有小血管但无中血管；

若高斯函数峰值小于第二设定阈值且方差小于第三设定阈值，则子图像第二设定灰度分布范围内既无小血管又无中血管；

若高斯函数峰值大于第二设定阈值且方差大于第三设定阈值，则子图像第二设定灰度分布范围内既有小血管又有中血管；

若高斯函数峰值小于第二设定阈值且方差大于第三设定阈值，则子图像第二设定灰度分布范围内无小血管但有中血管。

第三类：处于原图像边缘的子图像，编号为2、3、5、8、9、12、13、16、18、19。

1)同样，判断是否存在大血管。对位于原图边缘的子图像灰度直方图的第一设定灰度分布范围内像素点出现的频率进行统计，当统计值大于第四设定阈值时，则代表该子图像中含有大血管，反之，不含有大血管；但是由于这一类子图像分布在原图边缘，故而其中一定包含黑色背景，且这些黑色背景的像素点同样出现在直方图的第一设定灰度分布范围内。因此，这一类子图像用于判断大血管的阈值大于第二类子图像中用于判断大血管的阈值。

2)判断该子图像中是否包含小血管或中血管，与第二类子图像的判断标准相同。

若高斯函数峰值大于第五设定阈值且方差小于第六设定阈值，则子图像第二设定灰度分布范围内有小血管但无中血管；

若高斯函数峰值小于第五设定阈值且方差小于第六设定阈值，则子图像第二设定灰度分布范围内既无小血管又无中血管；

若高斯函数峰值大于第五设定阈值且方差大于第六设定阈值，则子图像第二设定灰度分布范围内既有小血管又有中血管；

若高斯函数峰值小于第五设定阈值且方差大于第六设定阈值，则子图像第二设定灰度分布范围内无小血管但有中血管。

步骤(4)：动态尺度分配：利用滤波器尺度与血管宽度之间的关系，动态选择不同尺度的滤波器对不同宽度血管进行增强；

所述步骤(4)的步骤为：

由于血管横截面的灰度分布呈高斯分布，故而将高斯曲线向横轴上投影，得到的投影长度就是血管的宽度，而正态曲线下[-3σ,+3σ]范围内的面积超过了99％，由此可以得出血管宽度与滤波器尺度之间的关系为：

σ＝d/6

其中，d表示血管宽度，σ表示滤波器的尺度。

a、本发明的仿真实验中，提取大血管所采用的滤波器尺度参数为2。

b、由各子图像的灰度直方图可以看出，中等宽度的血管与小血管在灰度分布范围上的界限并不明显。所以在二者灰度分布重叠的范围内，很可能会对中小血管的种类产生错误的判断。如若因此为其分配不恰当的尺度参数，将导致部分宽度的血管不能被检测出来。为避免这一情况的发生，本算法分别为提取中小血管分配了两个不同的尺度参数。这样，在整个算法中一共使用了五种尺度参数。以此在减小上述影响的同时，也使得更多宽度的血管得到增强。本发明仿真实验中，为提取中血管所采用的滤波器尺度参数为1和0.7，提取小血管所采用的滤波器尺度参数为0.5和0.2。

步骤(5)：多尺度匹配滤波；采用高斯匹配滤波模板和高斯一阶导数匹配滤波模板同时对增强后的视网膜图像进行匹配滤波，分别得到高斯匹配滤波器及高斯一阶导数匹配滤波器与视网膜图像的卷积响应；并进一步计算得到高斯一阶导数匹配滤波器与视网膜图像卷积响应的局部均值响应；

本发明采用两个匹配滤波模板，高斯型匹配滤波模板上坐标为(x,y)的一点，其权重表示如下：

其中，σ1表示高斯匹配滤波器的尺度，L代表与滤波器模板瞬时方向相平行的血管长度，m代表高斯滤波器模板中系数的均值，

式中Q表示构建的高斯滤波器模板中所包含的点数，N表示滤波器模板所在的邻域，表示高斯滤波器模板中的一个离散点。

在高斯型匹配滤波模板中减去m项，可以实现对视网膜图像背景的平滑处理。并且，由于视网膜血管在平面中的各个方向上均有分布，故而在匹配滤波过程中，滤波器模板需要不断的旋转以检测不同方向上的血管。

高斯一阶导数匹配滤波模板上坐标为(x,y)的一点，其权重表示如下：

其中，σ₂表示高斯一阶导数匹配滤波器的尺度，L代表与滤波器模板瞬时方向相平行的血管长度。

步骤(6)：阈值处理：将所有子图像利用其灰度分布特征，按照血管与背景的灰度对比度的高低，以及是否含具有台阶边缘的非血管结构进行分类，对各子图像进行动态的阈值处理；通过将整个图像的阈值矩阵与高斯匹配滤波器和视网膜图像的卷积响应逐点比较，每个像素点均与对应点的阈值进行比较，大于阈值的像素点被分类为血管，否则被分类为非血管；最后，将所有子图像的分割结果进行重新拼接，得到阈值处理后的视网膜血管网络二值图像；

阈值处理部分利用图像的灰度分布直方图，对每个子图像进行动态的阈值处理。将所有子图像利用其灰度分布特征按照对比度的高低及是否含有台阶边缘进行分类。

6-1)区分高对比度图像与低对比度图像：利用步骤(3)血管分类过程中得到的高斯曲线方差作为判断图像对比度高低的依据；

方差c小于设定值cth，表示子图像血管与背景的灰度对比度高，血管提取容易；

方差c大于设定值cth，表示子图像血管与背景的灰度对比度低，血管提取困难。

对低对比度子图像需进行如下处理：对第二设定灰度分布范围进行分段处理，以像素点出现频率最高的灰度值为界，分为高低两个灰度级范围。

6-2)区分是否存在台阶边缘：

对图像灰度直方图第三设定灰度分布范围内像素点出现的频率和s进行统计，并设立第七设定阈值，频率大于第七设定阈值时，代表子图像中含具有台阶边缘的非血管结构。

若判断出含具有台阶边缘的非血管结构，那么接下来要对其进行定位。进一步对局部均值响应的各点进行遍历判断，当中坐标为(x,y)的某一任意点响应幅度大于b_th时，认为(x,y)位置附近将会有具有台阶边缘的非血管结构像素点出现；小于b_th时则认为(x,y)位置附近将会有血管像素点出现。提取步骤(3)中所判定的三种不同规格血管时，所设置的b_th值也是不同的：

b_th＝1小血管

b_th＝1.1中血管

b_th＝1.2大血管

6-3)阈值选择

对各子图像的分割阈值由多阈值公式确定：

其中，T_k,j为子图像的阈值矩阵，h_i,j为常系数，为高斯滤波响应A的局部均值响应，k为图像血管与背景的灰度对比度对阈值个数的影响，j为有无台阶边缘对阈值个数的影响，K(c,c_th)与J(s,s_th)的取值具体如下：

其中，c是高斯拟合曲线的方差，当c小于阈值c_th时,表明对比度高，在整个灰度级范围内选择阈值，此时K的取值为1；反之，在对比度低的情况下，分为高低两个灰度级范围选择阈值，此时K的取值为2。s是第三设定灰度分布范围内像素点出现的频率和，如果s小于s_th,则不存在台阶边缘，不需对阈值针对血管结构和非血管结构进行区分，此时J的取值为1；当s大于等于s_th时,存在台阶边缘，血管结构和非血管结构要分别设置不同的阈值，此时J的取值为2。

针对以上分析，对各子图像的阈值处理主要分为下列四类情况：

第一类：对于对比度低且存在台阶边缘的子图像，采用h₁₁、h₁₂、h₂₁、h₂₂四个阈值，其中h₁₁针对高灰度级范围内的血管结构，h₁₂针对高灰度级范围内的非血管结构，h₂₁针对低灰度级范围内的血管结构，h₂₂针对低灰度级范围内的非血管结构；

第二类：对于对比度低且不存在台阶边缘的子图像，采用h₁₁、h₂₁两个阈值，其中h₁₁针对高灰度级范围内的血管结构，h₂₁针对低灰度级范围内的血管结构；

第三类：对于对比度高且存在台阶边缘的子图像，采用h₁₁、h₁₂两个阈值，其中h₁₁针对第二设定灰度分布范围内的血管结构，h₁₂针对第二设定灰度分布范围内的非血管结构；

第四类：对比度高且不存在台阶边缘，采用h₁₁一个阈值，h₁₁针对整个灰度级范围内的血管结构。

整个图像的阈值矩阵为最终可以通过下式得到分割后的血管网络：

其中，vess表示各像素点的最终取值。即将整幅图像的阈值矩阵与图像的高斯滤波响应A逐点比较，每个像素点均利用对应点的阈值判断一次，大于此点阈值的像素点被分类为血管，否则被分类为非血管。

最后，将各子图像的不同尺度分割结果利用逻辑“或”操作结合起来，再将20个子图像的分割结果进行重新拼接，得到阈值处理后的视网膜血管网络二值图像。

该动态阈值选择标准基于图像的灰度直方图，利用灰度直方图中体现的对比度和具有台阶边缘的非血管结构信息指导阈值的选择，与单阈值选择公式相比，既考虑了低对比度区域细小血管的提取，又考虑了对非血管结构的抑制，实验证明取得了更优的提取效果。

步骤(7)：后处理：对阈值处理后的视网膜血管网络二值图像进一步消除噪声，对阈值分割后的血管边缘进行平滑处理，消除图像中残留的视网膜边界，对不连续的细小血管进行断点连接，从而保留图像细节，得到分割精度高的视网膜血管网络图像。

首先，利用目标视网膜血管与噪点所在区域的几何特征差别，先删除视网膜血管网络二值图像中面积小于设定阈值P的连通域，再利用连通域外接矩形的长宽比，进一步的消除图像中残留的噪声；

其次，使用多尺度的高斯滤波器模版来对血管边缘进行平滑处理；

接下来，利用DRIVE数据库中提供的掩膜图像，找到视网膜边界像素点的位置，进而消除图像中残留的视网膜边界；

最后，采用闭运算对视网膜血管网络二值图像中不连续的细小血管进行断点连接。

本发明利用公共数据库DRIVE对提出的视网膜血管分割算法进行测试。此数据库中共有40张视网膜图像，是通过眼底照相机拍摄获得的，其中每一张图像的空间分辨率均为565×584。该数据库的建立者将数据库中的40张视网膜图像均等的分成了训练集和测试集，同时在这两个集合中均提供了一位专家手动分割的结果。但是对于测试集来说，它还提供了第二位专家的手动分割结果，并且，多数的算法均是采用测试集中的图像来进行仿真实验的。

将该数据库提供的第一位专家的手动分割结果，作为衡量本发明算法性能的标准(ground truth)。为了对不同视网膜血管分割算法的性能进行比较，定义了(1)TPR，(2)FPR，(3)ACC三个指标来对算法进行定量描述。定义如下：

其中ACC为视网膜血管分割平均准确率，TN为实验结果中，被正确分割的血管像素个数；BN为实验结果中，被正确分割的视场(FOV：图像中视网膜边界内的部分)中的背景像素个数。TPR定义为正确分类血管像素点的比例，Nvp为视网膜图像视场中，目标血管像素个数。FPR定义为视场内属于非血管但被分类为血管的像素点的比例，Nuvp为视网膜图像视场中，背景像素的个数。本发明对数据库DRIVE测试集中的全部20张图像进行了分割，取各性能指标的平均值作为本算法的最终结果，并与其他算法进行了性能比较，结果如表1所示。

表1是使用本发明算法对DRIVE数据库测试集全部20张图像进行处理，所得平均性能指标与其他算法的比较结果。图9(a)-图9(d)是本发明算法的20个提取结果之一。图9(a)为DRIVE数据库测试集中的第19号彩色视网膜血管图像，图9(b)为预处理后的图像，图9(c)为本发明算法提取的血管网络，图9(d)为ground truth血管网络图。对比图9(c)和图9(d)可以看出本发明算法能够将一些细小血管提取出来，并有效地抑制了非血管结构，具有较好的血管分割性能。

表1.DRIVE数据库测试集血管提取结果(with FOV)

如图2所示，基于动态尺度分配的视网膜血管提取系统，包括：

视网膜图像预处理模块：对彩色视网膜图像的绿色通道分量进行对比度增强；

图像分块模块：对预处理后的视网膜图像按照尺寸进行均等分块，分割成设定个数的子图像；

血管分类模块：根据血管宽度将血管分为大中小三类。并依照不同宽度血管的灰度分布特征，对子图像包含的血管种类进行判断；

动态尺度分配模块：利用滤波器尺度与血管宽度之间的关系，动态选择不同尺度的滤波器对不同宽度血管进行增强；

多尺度匹配滤波模块：采用高斯匹配滤波模板和高斯一阶导数匹配滤波模板同时对增强后的视网膜图像进行匹配滤波，分别得到高斯匹配滤波器与视网膜图像的卷积响应及高斯一阶导数匹配滤波器与视网膜图像的卷积响应；并进一步计算得到高斯一阶导数匹配滤波器与视网膜图像的局部均值响应；

阈值处理模块：利用高斯一阶导数匹配滤波器与视网膜图像卷积响应的局部均值响应对每个像素点的分割阈值进行调整，得到整个图像的阈值矩阵。将整个图像的阈值矩阵与高斯匹配滤波器和视网膜图像的卷积响应进行逐点比较，提取出血管结构并剔除非血管结构，将所有子图像的提取结果进行重新拼接，得到视网膜血管网络二值图像；

后处理模块：对阈值处理后的视网膜血管网络二值图像进一步消除噪声，对血管边缘进行平滑处理，消除图像中残留的视网膜边界，对不连续的细小血管进行断点连接，从而保留图像细节，得到分割精度高的视网膜血管网络图像。

上述虽然结合附图对本发明的具体实施方式进行了描述，但并非对本发明保护范围的限制，所属领域技术人员应该明白，在本发明的技术方案的基础上，本领域技术人员不需要付出创造性劳动即可做出的各种修改或变形仍在本发明的保护范围以内。