CN108961280B - 一种基于slic超像素分割的眼底视盘精细分割方法 - Google Patents

Abstract

本发明公开了一种基于SLIC超像素分割的眼底视盘精细分割方法，步骤为：对输入的眼底图像分别进行超像素分割、基于形态学处理的血管分割和R、G双通道颜色阈值分割，对颜色阈值分割后的连通域进行扩张处理后根据连通域的像素坐标选择相应的Toeplitz矩阵模板对眼底血管图像进行滤波，得到视盘中心位置。再提取视盘候选区域并去内部的血管，用二值化方法对视盘候选区域进行阈值分割并采用基于最小二乘的椭圆拟合方法确定视盘椭圆ROI区域，结合SLIC超像素分割结果保留具有一定重叠面积的超像素，实现视盘的精细分割。本发明实现视盘的自动定位和精细分割，可保留较好的视盘轮廓信息且耗时时间短，方便眼底图像的其他后续处理以及为眼科医生提供辅助诊断。

Description

一种基于SLIC超像素分割的眼底视盘精细分割方法

技术领域

本发明属于眼底图像视神经盘的分割方法技术领域，具体涉及一种基于SLIC超像素分割的眼底视盘精细分割方法。

背景技术

视神经盘即视盘，被认为是视网膜眼底图像的主要特征之一。视盘检测是设计用于视网膜解剖结构和病变的自动提取的许多算法中的一个关键预处理部分，也是大多数视网膜病变筛查系统的一个相关模块。视盘常作为其他眼底特征的标志物，如视盘与黄斑中央凹之间相当恒定的距离，可作为先验知识来帮助估计黄斑的位置。此外，视盘的形状、颜色或深度的变化是各种眼科病理学的指标，尤其是青光眼，因此视盘的维数常被用于测量某些视网膜病变引起的异常特征，例如青光眼和糖尿病视网膜病变等。此外，因为糖尿病视网膜病变的颜色与眼底渗出物相似，视盘作为虚警源经常会被被认为是一个或多个渗出候选区，因此正确的定位、分割和去除视盘可以提高渗出候选区的准确度，为渗出病变的检测打好基础。

视盘是正常眼底图像中最明亮的区域，它大致呈椭圆形。对视盘进行精确的定位与分割是十分困难的，特别是当眼底图像存在血管阻塞、边缘模糊和渗出性病变的干扰时。由于这些干扰通常与视神经盘有相似的亮度、颜色或对比度，有可能被误识别为视盘，因此所研究的各种方法应考虑不同图像之间的外观，大小和位置的变化，国内外学者基于此提出了多种视盘定位的方法。

Sinthanayothin等通过使用与视盘大小相同的窗函数对眼底图像进行滤波来识别相邻像素之间具有最高平均变化的区域，以此检测出视神经盘。Walter和Klein通过假设出现的明亮的病变区域大小与视盘大小相差很大，将视盘中心近似为眼底图像中最大的最亮连接对象的中心。Li和Chutatape通过亮度特征找到了基于模型的方法的视盘候选区域，他们选择灰度图像中具有最高1％灰度级的像素对其进行聚类，并丢弃小簇，并通过主成分分析(PCA)来创建视盘空间(OD模型)，再对候选区域中的每个像素运用具有不同比例的窗口进行PCA变换，最后视盘被检测为与其投影到盘空间上的欧式距离最小的区域。Lalonde通过金字塔分解采用多分辨率处理，对所有视盘候选区域计算简单置信度值后构建二进制边缘图，最后对边缘图和基于Hausdorff距离的模板匹配方法进行检测视盘。Frank terHaar通过搜索血管分支网络的所有路径，分析路径末端区域，并采用Hough变换来分割视盘。同时其他KNN、朴素贝叶斯和深度学习等方法也被用于视盘检测中。

检测视神经盘后，为了得到视神经盘更为细节的轮廓以便后续视神经盘具体情况的病变检测，需要对视神经盘进行精细分割。Walter等利用分水岭算法提取视神经盘轮廓。Lalonde等利用Canny算子检测出视盘边界，并将其与圆形模板进行匹配识别视盘区域。Lowell等提出一种基于全局椭圆模型和局部演化模型的方式进行视盘轮廓提取。此外，梯度矢量流(Gradient Vector Flow，GVF)、活动轮廓模型(Active Contour Model，ACM)、变分水平集模型等也被用于视盘轮廓的精细分割上。

然而，上述深度学习的网络训练依赖于庞大的样本数目；采用Hough变换或模板匹配等方法无法精细分割出视盘轮廓，对后续其他处理仍存在较大干扰；梯度矢量流、变分水平集模型等参数活动轮廓方法依赖于初始轮廓的设定，并且存在计算复杂度高和耗时长等缺点。因此，目前的视盘分割方法还存在轮廓分割不精细和算法耗时长等问题。

发明内容

本发明的目的在于：解决目前对眼底图像视盘精细分割时采用的方法依赖于庞大的样本数目或初始轮廓设定，计算复杂度高，耗时长且分割不精细的问题，提出一种基于SLIC超像素分割的眼底视盘精细分割方法，实现视盘的自动定位，采用最小二乘的椭圆拟合方法确定视盘椭圆ROI区域，结合SLIC超像素分割实现视盘的精细分割，可保留较好的视盘轮廓信息且耗时时间短，方便眼底图像的其他后续处理。

本发明采用的技术方案如下：

一种基于SLIC超像素分割的眼底视盘精细分割方法，方法步骤包括：

步骤1、对输入的眼底图像分别进行SLIC超像素分割、眼底血管图像分割和R、G双通道颜色阈值分割，对颜色阈值分割后的图像提取连通域，对连通域进行扩张处理再统计眼底图像垂直中心线两侧的像素数目；

步骤2、根据扩张后连通域图像垂直中心线两侧的像素数目选择相应的Toeplitz矩阵模板，并用模板对步骤1扩张后的连通区域和眼底血管图像做与操作后的图像进行滤波，得到眼底视盘定位位置坐标；

步骤3、根据眼底视盘定位位置坐标在输入的眼底图像原图中提取视盘候选区域，并去除视盘候选区域内部的血管；

步骤4、用二值化方法对步骤3得到的视盘候选区域进行阈值分割并采用基于最小二乘的椭圆拟合方法确定视盘椭圆ROI区域，基于步骤1SLIC超像素分割的结果，保留具有一定重叠面积的超像素，提取保留的超像素的轮廓，得到精细分割后的眼底视盘。

进一步，所述步骤1中SLIC超像素分割步骤具体为：

步骤1.1.11、对输入的眼底图像G首先将图像转化为五维特征向量V＝[l，a，b，x，y]，其中[l，a，b]为像素颜色，属于CIELAB颜色空间，[x，y]为像素位置；

步骤1.1.12、SLIC提出紧凑因子这种新的距离度量方法，即：

其中，k和i分别为两像素；D_s为CIELAB色彩空间值距离d_lab和图像平面内位置距离d_xy的加权和，表示两个像素间的距离；变量m度量超像素的紧凑性，m值越大，紧凑性就越高；若每幅图像像素的总数为N，预输出K个超像素，那么就有N/K个像素包含在每个超像素中，超像素的预期边长

且这些超像素在每个边长为S的网格中应有一个中心像素；

步骤1.1.13、基于紧凑因子的方法，将步骤1.1.11得到的眼底图像以网格为基本单位，在每个网格中选择一点作为超像素中心，计算其3*3邻域内像素的梯度，将梯度值最小的像素作为新的梯度中心；

步骤1.1.14、在每个梯度中心的2S×2S邻域内对属于该区域的像素进行搜索，并将所有像素归为与其临近的梯度中心，并对分割出的像素重新计算聚类中心，并计算新旧两区域中心的剩余误差；

步骤1.1.15、重复步骤1.1.12、1.1.13和1.1.14直到当剩余误差小于一定值或迭代次数超过一定值时，则超像素分割结束，得到超像素分割图像G_superpixel。

进一步，所述步骤1中基于形态学处理的眼底血管图像分割步骤具体为：

步骤1.1.21、首先提取输入的眼底图像G的G通道，并用中值滤波对其进行滤波，得到背景估计图像G_background，并用一定半径的圆盘结构元对背景估计图像G_background分别进行顶帽、底帽操作，得到G_top和G_black；

步骤1.1.22、将步骤1.1.21中的G_top和G_black以一定权重加在原图G上，并减去背景估计图像G_background，并进行阈值分割得到血管二值估计图G_BL；

步骤1.1.23、提取步骤1.1.22的血管二值估计图G_BL中连通域大于一定面积的连通区域，得到该眼底图像对应的眼底血管二值图像G_BV。

进一步，所述步骤1中对输入的眼底图像进行R、G双通道颜色阈值分割，对颜色阈值分割后的图像提取连通域并进行扩张处理步骤具体为：

步骤1.1.31、对输入眼底图像进行双通道，R通道和G通道阈值分割，得到二值图像G_RG；

步骤1.1.32、对步骤1.1.31中的二值图像G_RG进行连通域分析，提取每个连通域的最小矩形框，并对其进行向左、向右两个方向的扩张，得到图像G_expand；

步骤1.1.33、统计步骤1.1.32中的扩张连通域后图像G_expand中的图像垂直中心线两侧的像素数目。

进一步，所述步骤2的具体为：

步骤2.1、基于步骤1中眼底图像垂直中心线两侧的像素数目选择相应的Toeplitz矩阵模板，当统计的像素有超过半数位于图像右侧，则选择右Toeplitz矩阵模板Mask，如式(4)所示：

否则，选择左Toeplitz矩阵模板Mask，如式(5)所示：

步骤2.2、将步骤1中扩张后的连通区域和分割的眼底血管图像作与操作，得到候选区域的血管图像G_candidate：

G_candidate＝G_BV∩G_expand (6)；

步骤2.3、对候选区域的血管图像G_candidate进行滤波，得到滤波后的图像G_filter公式为：

步骤2.4、确定滤波后图像G_filter中灰度值最高的像素，该像素坐标即为眼底视盘定位位置坐标。

进一步，所述步骤3的具体为：

步骤3.1、以视盘定位位置坐标为中心，在输入的眼底图像G中框定一定大小矩形区域作为视盘候选区域；

步骤3.2、提取步骤3.1中视盘候选区域的R通道，记为图像G′_R；

步骤3.3、用半径递增的结构元进行交替膨胀腐蚀操作来去除视盘候选区域内部的血管，先选取一定半径的圆盘结构元B，对图像G′_R进行交替膨胀腐蚀操作：

σ^(B)＝δ^(B)(ε^(B)(G′_R)) (8)，

其中，δ^(B)表示以B为结构元进行膨胀操作；ε^(B)表示以B为结构元进行腐蚀操作；

步骤3.4、增加圆盘结构元B的半径，得到新的圆盘结构元B′，对步骤3.3处理后的图像再次进行交替膨胀腐蚀操作：

σ^(B′)＝δ^(B)(ε^(B′)(σ^(B))) (9)；

步骤3.5、继续增加圆盘结构元B′的半径，得到新的圆盘结构元B″，对步骤3.4处理后的图像再次进行交替膨胀腐蚀操作来去除视盘候选区域内部的血管：

σ^(B″)＝δ^(B″)(ε^(B″)(σ^(B))) (10)；

进一步，所述步骤4具体为：

步骤4.1、采用最大类间方差法对去除血管后的视盘候选区域进行阈值分割，得到二值图像G_Otsu，并用Canny算子提取二值图像G_Otsu的边缘，得到边缘图像G_edge；

步骤4.2、提取边缘图像G_edge中的轮廓坐标，对其进行最小二乘椭圆拟合；

步骤4.3、在眼底图像G绘制步骤4.2中椭圆的方程曲线，即为视盘椭圆ROI区域；

步骤4.4、在步骤1SLIC超像素分割后的图像中绘制步骤4.2中椭圆的方程曲线，并提取在该曲线内部的所有超像素集合S；

步骤4.5、计算步骤4.4中超像素集合S中每个超像素与视盘椭圆ROI区域的交集面积，若交集面积超过该超像素的一半面积，则保留该像素，否则舍弃该超像素；

步骤4.6、对步骤4.5保留的超像素进行轮廓提取，得到精细分割后的眼底视盘。

进一步，所述步骤4.2中最小二乘椭圆拟合的具体步骤为：

假设椭圆方程为：ax²+bxy+cy²+dx+ey＝1，则最小二乘椭圆拟合的最优化问题可以表示为：

min||Dα||²

s.t.α^TCα＝1 (11)，

其中α＝[a，b，c，d，e]；D表示轮廓坐标信息集合，维数为n×6，n为轮廓像素个数；矩阵C如下：

综上所述，由于采用了上述技术方案，本发明的有益效果是：

1、本发明中，解决目前对眼底图像视盘精细分割时采用的方法依赖于庞大的样本数目或初始轮廓设定，计算复杂度高，耗时长且分割不精细的问题，实现视盘的自动定位，采用最小二乘的椭圆拟合方法确定视盘椭圆ROI区域，结合SLIC超像素分割实现视盘的精细分割，可保留较好的视盘轮廓信息且耗时时间短，方便眼底图像的其他后续处理；

2、本发明中，采用最小二乘的椭圆拟合确定视盘ROI区域的方法，克服了常用的Hough变换的圆形拟合方法自适应差，在不同的图像可能需要不同的阈值的问题，且能更好的体现出视盘的形状，如果采用矩形框或者其他形状，结合后续的超像素分割之后，分出的视盘轮廓可能就不会很精确；

3、本发明中，对连通域进行扩张处理再统计眼底图像垂直中心线两侧的像素数目，根据扩张后连通域图像垂直中心线两侧的像素数目选择相应的Toeplitz矩阵模板，并用模板对扩张后的连通区域和眼底血管图像做与操作后的图像进行滤波，从而实现视盘的自动定位。连通域扩张处理解决了由于部分血管穿越视盘内部而导致双通道阈值分割无法分割出视盘中心区域的问题；相应的Toeplitz矩阵模板考虑了左右视盘的附近的血管走向分布，能更精确地实现视盘自动定位。

附图说明

为了更清楚地说明本发明实施例的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，应当理解，以下附图仅示出了本发明的某些实施例，因此不应被看作是对范围的限定，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他相关的附图。

图1为本发明方法流程图；

图2为本发明实施例1输入图像和血管分割结果；

图中(a)为输入眼底图像；(b)为血管分割结果；

图3为本发明实施例1视盘定位结果和视盘椭圆ROI；

图中(a)x符号标记处为视盘定位；(b)椭圆曲线标记处为x符号标记处对应的视盘椭圆ROI；

图4为本发明实施例1眼底图像超像素分割和视盘精细分割结果；

图中(a)为眼底图像超像素分割图像；(b)为视盘精细分割图像。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明，即所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。通常在此处附图中描述和示出的本发明实施例的组件可以以各种不同的配置来布置和设计。

因此，以下对在附图中提供的本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施例。基于本发明的实施例，本领域技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。

需要说明的是，术语“第一”和“第二”等之类的关系术语仅仅用来将一个实体或者操作与另一个实体或操作区分开来，而不一定要求或者暗示这些实体或操作之间存在任何这种实际的关系或者顺序。而且，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、物品或者设备所固有的要素。在没有更多限制的情况下，由语句“包括一个……”限定的要素，并不排除在包括所述要素的过程、方法、物品或者设备中还存在另外的相同要素。

一种基于SLIC超像素分割的眼底视盘精细分割方法，SLIC即simplelineariterative clustering：简单线性迭代聚类。本发明解决目前对眼底图像视盘精细分割时采用的方法依赖于庞大的样本数目或初始轮廓设定，计算复杂度高，耗时长且分割不精细的问题，实现视盘的自动定位，采用最小二乘的椭圆拟合方法确定视盘椭圆ROI区域，结合SLIC超像素分割实现视盘的精细分割，可保留较好的视盘轮廓信息且耗时时间短，方便眼底图像的其他后续处理。

方法步骤包括：

步骤l、对输入的眼底图像分别进行SLIC超像素分割、眼底血管图像分割和R、G双通道颜色阈值分割，对颜色阈值分割后的图像提取连通域，对连通域进行扩张处理再统计眼底图像垂直中心线两侧的像素数目；

步骤2、根据扩张后连通域图像垂直中心线两侧的像素数目选择相应的Toeplitz矩阵模板，并用模板对步骤1扩张后的连通区域和眼底血管图像做与操作后的图像进行滤波，得到眼底视盘定位位置坐标；

步骤3、根据眼底视盘定位位置坐标在输入的眼底图像原图中提取视盘候选区域，并去除视盘候选区域内部的血管；

步骤4、用二值化方法对步骤3得到的视盘候选区域进行阈值分割并采用基于最小二乘的椭圆拟合方法确定视盘椭圆ROI区域，基于步骤1SLIC超像素分割的结果，保留具有一定重叠面积的超像素，提取保留的超像素的轮廓，得到精细分割后的眼底视盘。

采用最小二乘的椭圆拟合确定视盘ROI区域的方法，克服了常用的Hough变换的圆形拟合方法自适应差，在不同的图像可能需要不同的阈值的问题，且尽可能好的体现出视盘的形状，如果采用矩形框或者其他形状，结合后续的超像素分割之后，分出的视盘轮廓可能就不会很精确。

本发明对连通域进行扩张处理再统计眼底图像垂直中心线两侧的像素数目，根据扩张后连通域图像垂直中心线两侧的像素数目选择相应的Toeplitz矩阵模板，并用模板对扩张后的连通区域和眼底血管图像做与操作后的图像进行滤波，从而实现视盘的自动定位。连通域扩张处理解决了由于部分血管穿越视盘内部而导致双通道阈值分割无法分割出视盘中心区域的问题；相应的Toeplitz矩阵模板考虑了左右视盘的附近的血管走向分布，能更精确地实现视盘自动定位。

进一步，所述步骤1中SLIC超像素分割步骤具体为：

步骤1.1.11、对输入的眼底图像G首先将图像转化为五维特征向量V＝[l，a，b，x，y]，其中[l，a，b]为像素颜色，属于CIELAB颜色空间，[x，y]为像素位置。

步骤1.1.12、由于颜色空间和距离空间的度量方法不同，SLIC提出紧凑因子这种新的距离度量方法，即：

其中，k和i分别为两像素；D_s为CIELAB色彩空间值距离d_lab和图像平面内位置距离d_xy的加权和，表示两个像素间的距离；变量m度量超像素的紧凑性，m值越大，紧凑性就越高；若每幅图像像素的总数为N，预输出K个超像素，那么就有N/K个像素包含在每个超像素中，超像素的预期边长

且这些超像素在每个边长为S的网格中应有一个中心像素。

步骤1.1.13、基于紧凑因子的方法，将步骤1.1.11得到的眼底图像以网格为基本单位，在每个网格中选择一点作为超像素中心，计算其3*3邻域内像素的梯度，将梯度值最小的像素作为新的梯度中心。

步骤1.1.14、在每个梯度中心的2S×2S邻域内对属于该区域的像素进行搜索，并将所有像素归为与其临近的梯度中心，并对分割出的像素重新计算聚类中心，并计算新旧两区域中心的剩余误差。

步骤1.1.15、重复步骤1.1.12、1.1.13和1.1.14直到当剩余误差小于一定值或迭代次数超过一定值时，则超像素分割结束，得到超像素分割图像G_superpixe。

进一步，所述步骤1中基于形态学处理的眼底血管图像分割步骤具体为：

步骤1.1.21、首先提取输入的眼底图像G的G通道，并用中值滤波对其进行滤波，得到背景估计图像G_background，并用一定半径的圆盘结构元对背景估计图像G_background分别进行顶帽、底帽操作，得到G_top和G_black。

步骤1.1.22、将步骤1.1.21中的G_top和G_black以一定权重加在原图G上，并减去背景估计图像G_background，并进行阈值分割得到血管二值估计图G_BL。

步骤1.1.23、提取步骤1.1.22的血管二值估计图G_BL中连通域大于一定面积的连通区域，得到该眼底图像对应的眼底血管二值图像G_BV。

进一步，所述步骤1中对输入的眼底图像进行R、G双通道颜色阈值分割，对颜色阈值分割后的图像提取连通域并进行扩张处理步骤具体为：

步骤1.1.31、对输入眼底图像进行双通道，R通道和G通道阈值分割，得到二值图像G_RG；

步骤1.1.32、对步骤1.1.31中的二值图像G_RG进行连通域分析，提取每个连通域的最小矩形框，并对其进行向左、向右两个方向的扩张，得到图像G_expand；

步骤1.1.33、统计步骤1.1.32中的扩张连通域后图像G_expand中的图像垂直中心线两侧的像素数目。

进一步，所述步骤2的具体为：

步骤2.1、基于步骤1中眼底图像垂直中心线两侧的像素数目选择相应的Toeplitz矩阵模板，当统计的像素有超过半数位于图像右侧，则选择右Toeplitz矩阵模板Mask，如式(4)所示：

否则，选择左Toeplitz矩阵模板Mask，如式(5)所示：

步骤2.2、将步骤1中扩张后的连通区域和分割的眼底血管图像作与操作，得到候选区域的血管图像G_candidate：

G_candidate＝G_BV∩G_expand (6)。

步骤2.3、对候选区域的血管图像G_candidate进行滤波，得到滤波后的图像G_filter公式为：

步骤2.4、确定滤波后图像G_filter中灰度值最高的像素，该像素坐标即为眼底视盘定位位置坐标。

进一步，所述步骤3的具体为：

步骤3.1、以视盘定位位置坐标为中心，在输入的眼底图像G中框定一定大小矩形区域作为视盘候选区域。

步骤3.2、提取步骤3.1中视盘候选区域的R通道，记为图像G′_R。

步骤3.3、用半径递增的结构元进行交替膨胀腐蚀操作来去除视盘候选区域内部的血管，也可以用其他方法代替，只要能将视盘候选区域内部血管去掉的方法都可以，这里用半径递增的结构元进行交替膨胀腐蚀操作进行说明，先选取一定半径的圆盘结构元B，对图像G′_R进行交替膨胀腐蚀操作：

σ^(B)＝δ^(B)(ε^(B)(G′_R)) (8)，

其中，δ^(B)表示以B为结构元进行膨胀操作；ε^(B)表示以B为结构元进行腐蚀操作。

步骤3.4、增加圆盘结构元B的半径，得到新的圆盘结构元B′，对步骤3.3处理后的图像再次进行交替膨胀腐蚀操作：

σ^(B′)＝δ^(B′)(ε^(B′)(σ^(B))) (9)。

步骤3.5、继续增加圆盘结构元B′的半径，得到新的圆盘结构元B″，对步骤3.4处理后的图像再次进行交替膨胀腐蚀操作来去除视盘候选区域内部的血管：

σ^(B″)＝δ^(B″)(ε^(B″)(σ^(B′))) (10)。

进一步，步骤4中可采用最大类间方差法对步骤3得到的视盘候选区域进行阈值分割，最大类间方差法可以自动寻找到一个阈值来将灰度图像二值化，但也可以人工设定阈值，或者采用其他方法进行二值化，这里采用最大类间方差法详细说明步骤4具体为：

步骤4.1、采用最大类间方差法对步骤3.5得到的去除血管后的视盘候选区域进行阈值分割，得到二值图像G_Otsu，并用Canny算子提取二值图像G_Otsu的边缘，得到边缘图像G_edge。

步骤4.2、提取边缘图像G_edge中的轮廓坐标，对其进行最小二乘椭圆拟合。

所述步骤4.2中最小二乘椭圆拟合的具体步骤为：

假设椭圆方程为：ax²+bxy+cy²+dx+ey＝1，则最小二乘椭圆拟合的最优化问题可以表示为：

min||Dα||²

s.t.α^TCα＝1 (11)，

其中α＝[a，b，c，d，e]；D表示轮廓坐标信息集合，维数为n×6，n为轮廓像素个数；矩阵C如下：

步骤4.3、在眼底图像G绘制步骤4.2中椭圆的方程曲线，即为视盘椭圆ROI区域。

步骤4.4、在步骤1SLIC超像素分割后的图像中绘制步骤4.2中椭圆的方程曲线，并提取在该曲线内部的所有超像素集合S。

步骤4.5、计算步骤4.4中超像素集合S中每个超像素与视盘椭圆ROI区域的交集面积，若交集面积超过该超像素的一半面积，则保留该像素，否则舍弃该超像素。

步骤4.6、对步骤4.5保留的超像素进行轮廓提取，得到精细分割后的眼底视盘。

以下结合实施例对本发明的特征和性能作进一步的详细描述。

实施例1

本发明较佳实施例提供的一种基于SLIC超像素分割的眼底视盘精细分割方法，流程图如图1所示，方法步骤包括：

步骤1、对输入的眼底图像分别进行SLIC超像素分割、眼底血管图像分割和R、G双通道颜色阈值分割，对颜色阈值分割后的图像提取连通域，对连通域进行扩张处理再统计眼底图像垂直中心线两侧的像素数目。

所述步骤1中SLIC超像素分割步骤具体为：

步骤1.1.11、对输入的眼底图像G首先将图像转化为五维特征向量V＝[l，a，b，x，y]，其中[l，a，b]为像素颜色，属于CIELAB颜色空间，[x，y]为像素位置。

步骤1.1.12、由于颜色空间和距离空间的度量方法不同，SLIC提出紧凑因子这种新的距离度量方法，即：

其中，k和i分别为两像素；D_s为CIELAB色彩空间值距离d_lab和图像平面内位置距离d_xy的加权和，表示两个像素间的距离；变量m度量超像素的紧凑性，m值越大，紧凑性就越高；若每幅图像像素的总数为N，预输出K个超像素，那么就有N/K个像素包含在每个超像素中，超像素的预期边长

且这些超像素在每个边长为S的网格中应有一个中心像素。本实施例中，输入图像大小为1552×1928，预输出220个超像素，且紧凑系数m＝15。

步骤1.1.13、基于紧凑因子的方法，将步骤1.1.11得到的眼底图像以网格为基本单位，在每个网格中选择一点作为超像素中心，计算其3*3邻域内像素的梯度，将梯度值最小的像素作为新的梯度中心。

步骤1.1.14、在每个梯度中心的2S×2S邻域内对属于该区域的像素进行搜索，并将所有像素归为与其临近的梯度中心，并对分割出的像素重新计算聚类中心，并计算新旧两区域中心的剩余误差。

步骤1.1.15、重复步骤1.1.12、1.1.13和1.1.14直到当剩余误差小于一定值或迭代次数超过一定值时，则超像素分割结束，得到超像素分割图像G_superpixel。

所述步骤1中基于形态学处理的眼底血管图像分割步骤具体为：

步骤1.1.21、首先提取输入的眼底图像G的G通道，并用中值滤波对其进行滤波，得到背景估计图像G_background，并用一定半径的圆盘结构元对背景估计图像G_background分别进行顶帽、底帽操作，得到G_top和G_black。本实施例中，采用40×40大小的中值滤波器，半径为20的圆盘结构元。

步骤1.1.22、将步骤1.1.21中的G_top和G_black以一定权重加在原图G上，并减去背景估计图像G_background，并进行阈值分割得到血管二值估计图G_BL。本实施例中得到G_BL的公式为：G_BL＝l.3*G_top-1.2*G_black+G-G_background。

步骤1.1.23、提取步骤1.1.22的血管二值估计图G_BL中连通域大于一定面积的连通区域，得到该眼底图像对应的眼底血管二值图像G_BV。本示例中，连通区域面积阈值取为400个像素值。图2(a)为输入眼底图像，(b)为血管分割的结果。

进一步，所述步骤1中对输入的眼底图像进行R、G双通道颜色阈值分割，对颜色阈值分割后的图像提取连通域并进行扩张处理步骤具体为：

步骤1.1.31、对输入眼底图像进行双通道，R通道和G通道阈值分割，得到二值图像G_RG。RG通道各像素值范围为0-255，本实施例中，R通道阈值为250，G通道阈值为160。

步骤1.1.32、统计步骤1.1.31中的二值图像G_RG中的图像垂直中心线两侧的像素数目。

步骤1.1.33、对步骤1.1.31中的二值图像G_RG进行连通域分析，提取每个连通域的最小矩形框，并对其进行向左、向右两个方向的扩张一倍最小矩形框大小的距离，得到图像G_expand。

步骤2、根据扩张后连通域图像垂直中心线两侧的像素数目选择相应的Toeplitz矩阵模板，并用模板对步骤1扩张后的连通区域和眼底血管图像做与操作后的图像进行滤波，得到眼底视盘定位位置坐标。

所述步骤2的具体为：

步骤2.1、基于步骤1中眼底图像垂直中心线两侧的像素数目选择相应的Toeplitz矩阵模板，当统计的像素有超过半数位于图像右侧，则选择右Toeplitz矩阵模板Mask，如式(4)所示：

否则，选择左Toeplitz矩阵模板Mask，如式(5)所示：

步骤2.2、将步骤1中扩张后的连通区域和分割的眼底血管图像作与操作，得到候选区域的血管图像G_candidate：

G_candidate＝G_BV∩G_expand (6)。

步骤2.3、对候选区域的血管图像G_candidate进行滤波，得到滤波后的图像G_filter公式为：

步骤2.4、确定滤波后图像G_filter中灰度值最高的像素，该像素坐标即为眼底视盘定位位置坐标。

步骤3、根据眼底视盘定位位置坐标在输入的眼底图像原图中提取视盘候选区域，并用半径递增的结构元进行交替膨胀腐蚀操作来去除视盘候选区域内部的血管。

所述步骤3的具体为：

步骤3.1、以视盘定位位置坐标为中心，在输入的眼底图像G中框定一定大小矩形区域作为视盘候选区域。本实施例中矩形框大小为400×400。

步骤3.2、提取步骤3.1中视盘候选区域的R通道，记为图像G′_R。

步骤3.3、选取一定半径的圆盘结构元B，对图像G′_R进行交替膨胀腐蚀操作：

σ^(B)＝δ^(B)(ε^(B)(G′_R)) (8)，

其中，δ^(B)表示以B为结构元进行膨胀操作；ε^(B)表示以B为结构元进行腐蚀操作。

步骤3.4、增加圆盘结构元B的半径，得到新的圆盘结构元B′，对步骤3.3处理后的图像再次进行交替膨胀腐蚀操作：

σ^(B′)＝δ^(B′)(ε^(B′)(σ^(B))) (9)。

步骤3.5、继续增加圆盘结构元B′的半径，得到新的圆盘结构元B″，对步骤3.4处理后的图像再次进行交替膨胀腐蚀操作来去除视盘候选区域内部的血管：

σ^(B″)＝δ^(B″)(ε^(B″)(σ^(B′))) (10)。

步骤4、用最大类间方差法对步骤3得到的视盘候选区域进行阈值分割并采用椭圆拟合方法确定视盘椭圆ROI区域，基于步骤1SLIC超像素分割的结果，保留具有一定重叠面积的超像素，提取保留的超像素的轮廓，得到精细分割后的眼底视盘。图3为视盘定位结果和视盘椭圆ROI，(a)中x符号标记处为视盘定位，(b)中椭圆曲线标记处为x对应的视盘椭圆ROI。

所述步骤4的具体为：

步骤4.1、采用最大类间方差法对步骤3得到的视盘候选区域进行阈值分割，得到二值图像G_Otsu，并用Canny算子提取二值图像G_Otsu的边缘，得到边缘图像G_edge。

步骤4.2、提取边缘图像G_edge中的轮廓坐标，对其进行最小二乘椭圆拟合。

所述步骤4.2中最小二乘椭圆拟合的具体步骤为：

假设椭圆方程为：ax²+bxy+cy²+dx+ey＝1，则最小二乘椭圆拟合的最优化问题可以表示为：

min||Dα||²

s.t.α^TCα＝1 (11)，

其中α＝[a，b，c，d，e]；D表示轮廓坐标信息集合，维数为n×6，n为轮廓像素个数；矩阵C如下：

步骤4.3、在眼底图像G绘制步骤4.2中椭圆的方程曲线，即为视盘椭圆ROI区域。

步骤4.4、在步骤1SLIC超像素分割后的图像中绘制步骤4.2中椭圆的方程曲线，并提取在该曲线内部的所有超像素集合S。

步骤4.5、计算步骤4.4中超像素集合S中每个超像素与视盘椭圆ROI区域的交集面积，若交集面积超过该超像素的一半面积，则保留该像素，否则舍弃该超像素。

步骤4.6、对步骤4.5保留的超像素进行轮廓提取，得到精细分割后的眼底视盘。

图4为眼底图像超像素分割和视盘精细分割结果，(a)为眼底图像超像素分割示意图，(b)为视盘精细分割示意图。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (8)

1.一种基于SLIC超像素分割的眼底视盘精细分割方法，其特征在于：方法步骤包括：

步骤1、对输入的眼底图像分别进行SLIC超像素分割、基于形态学处理的眼底血管图像分割和R、G双通道颜色阈值分割，对颜色阈值分割后的图像提取连通域，对连通域进行扩张处理再统计眼底图像垂直中心线两侧的像素数目；

步骤2、根据扩张后连通域图像垂直中心线两侧的像素数目选择相应的Toeplitz矩阵模板，并用模板对步骤1扩张后的连通区域和眼底血管图像做与操作后的图像进行滤波，得到眼底视盘定位位置坐标；

步骤3、根据眼底视盘定位位置坐标在输入的眼底图像原图中提取视盘候选区域，并去除视盘候选区域内部的血管；

步骤4、用二值化方法对步骤3得到的视盘候选区域进行阈值分割并采用基于最小二乘的椭圆拟合方法确定视盘椭圆ROI区域，基于步骤1的SLIC超像素分割的结果，保留具有一定重叠面积的超像素，提取保留的超像素的轮廓，得到精细分割后的眼底视盘。

2.根据权利要求1所述的一种基于SLIC超像素分割的眼底视盘精细分割方法，其特征在于：步骤1中所述SLIC超像素分割步骤具体为：

步骤1.1.11、对输入的眼底图像G首先将图像转化为五维特征向量V＝[l,a,b,x,y]，其中[l,a,b]为像素颜色，属于CIELAB颜色空间，[x,y]为像素位置；

步骤1.1.12、SLIC提出紧凑因子这种新的距离度量方法，即：

其中，k和i分别为两像素；D_s为CIELAB色彩空间值距离d_lab和图像平面内位置距离d_xy的加权和，表示两个像素间的距离；变量m度量超像素的紧凑性，m值越大，紧凑性就越高；若每幅图像像素的总数为N，预输出K个超像素，那么就有N/K个像素包含在每个超像素中，超像素的预期边长

且这些超像素在每个边长为S的网格中应有一个中心像素；

步骤1.1.13、基于紧凑因子的方法，将步骤1.1.11得到的眼底图像以网格为基本单位，在每个网格中选择一点作为超像素中心，计算其3*3邻域内像素的梯度，将梯度值最小的像素作为新的梯度中心；

步骤1.1.14、在每个梯度中心的2S×2S邻域内对属于该区域的像素进行搜索，并将所有像素归为与其临近的梯度中心，并对分割出的像素重新计算聚类中心，并计算新旧两区域中心的剩余误差；

步骤1.1.15、重复步骤1.1.12、1.1.13和1.1.14直到当剩余误差小于一定值或迭代次数超过一定值时，则超像素分割结束，得到超像素分割图像G_superpixel。

3.根据权利要求1所述的一种基于SLIC超像素分割的眼底视盘精细分割方法，其特征在于：步骤1中所述基于形态学处理的眼底血管图像分割步骤具体为：

步骤1.1.21、首先提取输入的眼底图像G的G通道，并用中值滤波对其进行滤波，得到背景估计图像G_background，并用一定半径的圆盘结构元对背景估计图像G_background分别进行顶帽、底帽操作，得到G_top和G_black；

步骤1.1.22、将步骤1.1.21中的G_top和G_black以一定权重加在原图G上，并减去背景估计图像G_background,并进行阈值分割得到血管二值估计图G_BL；

步骤1.1.23、提取步骤1.1.22的血管二值估计图G_BL中连通域大于一定面积的连通区域，得到该眼底图像对应的眼底血管二值图像G_BV。

4.根据权利要求1所述的一种基于SLIC超像素分割的眼底视盘精细分割方法，其特征在于：所述步骤1中对输入的眼底图像进行R、G双通道颜色阈值分割，对颜色阈值分割后的图像提取连通域并进行扩张处理步骤具体为：

步骤1.1.31、对输入眼底图像进行双通道，R通道和G通道阈值分割，得到二值图像G_RG；

步骤1.1.32、对步骤1.1.31中的二值图像G_RG进行连通域分析，提取每个连通域的最小矩形框，并对其进行向左、向右两个方向的扩张，得到图像G_expand；

步骤1.1.33、统计步骤1.1.32中的扩张连通域后图像G_expand中的图像垂直中心线两侧的像素数目。

5.根据权利要求1所述的一种基于SLIC超像素分割的眼底视盘精细分割方法，其特征在于：所述步骤2中的具体步骤为：

步骤2.1、基于步骤1中眼底图像垂直中心线两侧的像素数目选择相应的Toeplitz矩阵模板，当统计的像素有超过半数位于图像右侧，则选择右Toeplitz矩阵模板Mask，如式(4)所示：

否则，选择左Toeplitz矩阵模板Mask，如式(5)所示：

步骤2.2、将步骤1中扩张后的连通区域和分割的眼底血管图像作与操作，得到候选区域的血管图像G_candidate：

G_candidate＝G_BV∩G_expand (6)；

步骤2.3、对候选区域的血管图像G_candidate进行滤波，得到滤波后的图像G_filter公式为：

步骤2.4、确定滤波后图像G_filter中灰度值最高的像素，该像素坐标即为眼底视盘定位位置坐标。

6.根据权利要求1所述的一种基于SLIC超像素分割的眼底视盘精细分割方法，其特征在于：所述步骤3中的具体步骤为：

步骤3.1、以视盘定位位置坐标为中心，在输入的眼底图像G中框定一定大小矩形区域作为视盘候选区域；

步骤3.2、提取步骤3.1中视盘候选区域的R通道，记为图像G_R′；

步骤3.3、用半径递增的结构元进行交替膨胀腐蚀操作来去除视盘候选区域内部的血管，先选取一定半径的圆盘结构元B，对图像G_R′ 进行交替膨胀腐蚀操作：

σ^(B)＝δ^(B)(ε^(B)(G_R′)) (8)，

其中，δ^(B)表示以B为结构元进行膨胀操作；ε^(B)表示以B为结构元进行腐蚀操作；

步骤3.4、增加圆盘结构元B的半径，得到新的圆盘结构元B′，对步骤3.3处理后的图像再次进行交替膨胀腐蚀操作：

σ^(B′)＝δ^(B′)(ε^(B′)(σ^(B))) (9)；

步骤3.5、继续增加圆盘结构元B′的半径，得到新的圆盘结构元B″，对步骤3.4处理后的图像再次进行交替膨胀腐蚀操作来去除视盘候选区域内部的血管：

σ^(B″)＝δ^(B″)(ε^(B″)(σ^(B′))) (10)。

7.根据权利要求1所述的一种基于SLIC超像素分割的眼底视盘精细分割方法，其特征在于：所述步骤4具体为：

步骤4.1、采用最大类间方差法对去除血管后的视盘候选区域进行阈值分割，得到二值图像G_Otsu，并用Canny算子提取二值图像G_Otsu的边缘，得到边缘图像G_edge；

步骤4.2、提取边缘图像G_edge中的轮廓坐标，对其进行最小二乘椭圆拟合；

步骤4.3、在眼底图像G绘制步骤4.2中椭圆的方程曲线，即为视盘椭圆ROI区域；

步骤4.4、在步骤1SLIC超像素分割后的图像中绘制步骤4.2中椭圆的方程曲线，并提取在该曲线内部的所有超像素集合S；

步骤4.5、计算步骤4.4中超像素集合S中每个超像素与视盘椭圆ROI区域的交集面积，若交集面积超过该超像素的一半面积，则保留该像素，否则舍弃该超像素；

步骤4.6、对步骤4.5保留的超像素进行轮廓提取，得到精细分割后的眼底视盘。

8.根据权利要求7所述的一种基于SLIC超像素分割的眼底视盘精细分割方法，其特征在于：所述步骤4.2中最小二乘椭圆拟合的具体步骤为：

假设椭圆方程为：ax²+bxy+cy²+dx+ey＝1，则最小二乘椭圆拟合的最优化问题可以表示为：

min||Dα||²

s.t.α^TCα＝1 (11)，

其中α＝[a,b,c,d,e]；D表示轮廓坐标信息集合，维数为n×6，n为轮廓像素个数；矩阵C如下：

Images