CN108986106B - 面向青光眼的视网膜血管自动分割方法 - Google Patents

Abstract

本发明提出一种面向青光眼临床诊断的视网膜血管自动分割方法，通过融合匹配滤波器、神经网络、多尺度线检测、尺度空间分析和形态学五种依赖于不同图像处理技术的模型，消除了视盘和渗出物等明亮区域的影响。同时，本发明不需要海量级别的数据建立视网膜血管分割模型，极大地降低了需要处理的数据量和复杂度，易于实现，能有效提高视网膜血管分割的效率。本发明还在多模态融合结果的基础上利用区域生长法和梯度信息对背景和血管区域进行迭代生长，分割结果具有较好的连续性和平滑性，能保留更多的视网膜血管细节和更完整视的网膜血管网络，从而有效协助眼科医生诊断疾病，减轻眼科医生的负担。

Description

面向青光眼的视网膜血管自动分割方法

技术领域

本发明涉及数字医学图像处理与分析和智能医疗领域，具体涉及一种面向青光眼临床诊断的视网膜血管自动分割方法。

背景技术

青光眼是全球首要不可逆性致盲眼病，被称为“无声的视力窃贼”。预计到2020年，全球青光眼患病人数将会上升到7960万人。目前我国40岁以上人群开角型青光眼的患病率达到2.6％，占青光眼患者总人数的2/3，致盲率为15％～30％，远高于经济发达国家8％的平均水平。早发现、早诊断、早治疗对抑制青光眼病情的发展十分重要。彩色眼底图像能直接观测到视网膜血管病变和其他诸如渗出物、微动脉瘤等病灶，被眼科医生广泛应用于青光眼的临床诊断。视网膜血管是彩色眼底图像中可见的最主要解剖结构，其结构特征变化能直接反映青光眼、高血压等血管性相关疾病对血管网络形态结构的影响。视网膜血管的分割对青光眼的筛查、诊断和治疗等至关重要，是青光眼诊治的必要步骤。临床上诊断青光眼由有经验的眼科医生完成，首先手动分割患者的视网膜血管，然后再测量血管管径、分叉角度等所需相关参数。其中，手动分割视网膜血管的工作繁琐，需花费大量的时间和精力。

视网膜血管自动分割是一项极具挑战性的工作，因为虽然视网膜血管和背景有一定区别，但其亮度会随着血管的延伸而逐渐变化，尤其是血管末梢和背景的对比度低，大大增加了视网膜血管完全分割的难度。目前已有的视网膜血管自动分割方法大致分为非监督分割法和监督学习分割法两大类。非监督分割法无需真值图作为金标准对模型进行训练，主要包括局部自适应阈值法、局部滤波法和血管追踪法等。监督学习分割法将真值图数据作为训练样本对模型进行训练，逐步优化模型参数后分割视网膜血管；主要包括基于神经网络的像素分类法和基于脊波的血管分割法等。专利号为201710418711.4的发明专利《糖尿病人视网膜血管图像分割方法》，包括：按照预设的二值化阈值对经过预处理的眼底图像进行二值化处理，并提取二值化处理后的眼底图像中的中心线和边缘，得到血管树；对所述的血管树分叉处做断开处理得到血管段，并对各个血管段进行线分割得到血管，得到原始血管集；确定误分割血管，并从原始血管集中去除得到全局血管集。这些视网膜血管分割方法虽然取得了一定的研究成果，但仍存在一些关键问题亟需解决，主要存在分割出的微细视网膜血管易离散、连续性不佳以及难以有效地提取完整血管网络等局限性，应用范围有限。

近来，深度学习在语音识别、目标检测、图像识别和遥感图像分类等应用领域获得了很大进展，成为当前的研究热点之一。专利号为201610844032.9的发明专利《基于深度学习的眼底图像视网膜血管分割方法及系统》，包括：对训练集进行数据扩增，并对图像进行增强，用训练集训练卷积神经网络，先使用卷积神经网络分割模型对图像进行分割得到一个分割结果，用卷积神经网络的特征训练随机森林分类器，从卷积神经网络模型中抽取最后一层卷积层输出，并作为随机森林分类器的输入进行像素分类，得到另外一个分割结果，对两个分割结果进行融合得到最终的分割图像。与传统方法相比，该方法虽然能够获得更好的分割准确率和鲁棒性；但是，深度学习是一种数据驱动型模型，需要海量级别的数据作保证，这会严重影响视网膜血管分割的效率，限制了其在临床实践的应用推广。

说明书中提到Zhu、Chaudhuri、Mendonca、Staal、Soares、Zhang、Fraz、Zhao、Franklin和Wang多种方法与本发明进行对比，各方法对应的公开文献分别如下：

Zhu C Z,Xiang Y,Zou B J,et al.Retinal vessel segmentation in fundusimages using CART and AdaBoost[J].Journal of Computer-Aided Design&ComputerGraphics,2014,26(3):445-451.

Chaudhuri S,Chatterjee S,Katz N,et al.Detection of blood vessels inretinal images using two-dimensional matched filters[J].IEEE Transactions onMedical Imaging,1989,8(3):263-269.

Mendonca A M,Campilho A.Segmentation of retinal blood vessels bycombining the detection of centerlines and morphological reconstruction[J].IEEE Transactions on Medical Imaging,2006,25(9):1200-1213.

Staal J,Abràmoff M D,Niemeijer M,et al.Ridge-based vesselsegmentation in color images of the retina[J].IEEE Transactions on MedicalImaging,2004,23(4):501-509.

Soares J V,Leandro J J G,Cesar R M,et al.Retinal vessel segmentationusing the 2-D Gabor wavelet and supervised classification[J].IEEETransactions on Medical Imaging,2006,25(9):1214-1222.

Zhang B,Zhang L,Zhang L,et al.Retinal vessel extraction by matchedfilter with first-order derivative of Gaussian[J].Computers in Biology andMedicine,2010,40(4):438-445.

Fraz M M,Barman S A,Remagnino P,et al.An approach to localize theretinal blood vessels using bit planes and centerline detection[J].ComputerMethods and Programs in Biomedicine,2012,108(2):600-616.

Zhao Y Q,Wang X H,Wang X F,et al.Retinal vessels segmentation basedon level set and region growing[J].Pattern Recognition,2014,47(7):2437-2446.

Franklin S W,Rajan S E.Computerized screening of diabetic retinopathyemploying blood vessel segmentation in retinal images[J].Biocybernetics andBiomedical Engineering,2014,34(2):117-124.

Wang S,Yin Y,Cao G,et al.Hierarchical retinal blood vesselsegmentation based on feature and ensemble learning[J].Neurocomputing,2015,149(B):708-717.

发明内容

本发明要解决的技术问题是提出一种通过多模型融合实现分割的、面向青光眼临床诊断的视网膜血管自动分割方法。

为解决上述技术问题，本发明提出面向青光眼临床诊断的视网膜血管自动分割方法，包括以下步骤：

步骤1：眼底图像预处理：预处理待分割眼底图像；

步骤2：初步分割：分别构建匹配滤波器模型、神经网络模型、多尺度线检测模型、尺度空间分析模型和形态学模型；利用上述构建的模型，分别对步骤1所得待分割眼底图像进行初步分割，获得对应的初步分割结果；将五个初步分割结果的均值作为初步分割输出；

步骤3：多模型融合：利用掩膜分离步骤1所得待分割眼底图像中渗出物和视盘区域，并利用步骤2所得的形态学模型分割结果替换掩膜中白色区域，再与步骤2所得的初步分割输出融合，生成组合结果；

步骤4：精分割：利用Otsu法阈值分割步骤3所得组合结果，并根据血管的连通特性进行区域迭代生长后，获取视网膜血管分割的最终结果。

作为本发明面向青光眼临床诊断的视网膜血管自动分割方法的改进：

所述步骤1中眼底图像预处理针对匹配滤波器模型、神经网络模型、多尺度线检测模型、尺度空间分析模型和形态学模型采取不同预处理步骤：

1.1、匹配滤波器模型、尺度空间分析模型以及形态学分割模型的预处理：

将待分割眼底图像分解为红色、绿色和蓝色三个分量图像，提取待分割眼底图像的绿色分量图像；

1.2、神经网络分割模型的预处理：

按照步骤1.1提取待分割眼底图像的绿色分量图像，并对其进行以下处理：

1.2.1、利用圆盘状结构元素(直径为3像素)对绿色分量图像进行形态学“开”运算，并对所得图像进行均值滤波(窗口大小为69×69像素)，得到背景图像；

1.2.2、利用按照步骤1.1所提取的绿色分量图像减去步骤1.2.1所得的背景图像，并将像素灰度值调整至0-1范围内实现背景均匀化，获得背景均匀化后的图像；

1.2.3、首先将步骤1.2.2所得背景均匀化后的图像中血管进行增强，之后对该图像进行“补”运算，并使用圆盘状结构元素(半径为8像素)对“补”运算后的图像进行top-hat(顶帽)变换，增强该图像的暗区域以及视网膜血管区域，同时去除明亮区域(如视盘区域)，获得血管增强后的图像；

1.3、多尺度线检测分割模型的预处理：

首先按照步骤1.1中提取眼底图像的绿色分量图像，并对该绿色分量图像进行“取反”操作，之后对所得图像进行背景均匀化处理；之后将处理后的图像进行均值滤波(窗口大小为69×69像素)，并将均值滤波后的局部平均灰度替换明亮区域。

作为本发明面向青光眼临床诊断的视网膜血管自动分割方法的进一步改进：

所述步骤2中构建匹配滤波器分割模型步骤如下：

所述匹配滤波器分割模型为二维匹配滤波器，其高斯曲线的表达式为：

其中，K(x,y)称为高斯核函数；(x,y)表示步骤1.1所得绿色分量图像中像素的坐标，满足|x|≤λσ，|y|≤L/2；σ为高斯核沿x轴坐标中心的离散度；L为沿y轴被截断视网膜血管的长度；λ为常数；

所述高斯核函数为旋转矩阵(11×11)，在0°～180°范围内，每隔15°就选取一个高斯核函数，共选取了12个不同的高斯核旋转矩阵；

将步骤1.1所得绿色分量图像通过匹配滤波器分割模型进行分割的步骤如下：

将步骤1.1高斯核旋转矩阵中每个像素分别与12个不同的高斯核旋转矩阵进行卷积，选取其中最大的卷积值作为匹配滤波后的像素值，得到视网膜血管概率图；之后对视网膜血管概率图进行阈值分割(全局阈值t_mf设置为0.33)，得到含有背景区域及非背景区域的二值化图像；

将所处区域面积大于950像素，且该区域为非背景区域的像素标记为视网膜血管像素，合并所有视网膜血管像素，获得匹配滤波器分割模型的初步分割结果。

作为本发明面向青光眼临床诊断的视网膜血管自动分割方法的进一步改进：

所述步骤2中构建神经网络分割模型的步骤如下：

所述馈反向传播神经网络为多层的前馈反向传播神经网络分割模型，包含1个输入层(7个节点)、3个隐层(15个节点)和1个输出层(1个节点)；隐层的传输函数采用线性函数，输出层的传输函数logsig(z)采用对数S型函数：

其中，z表示隐层的输出值；

所述前馈反向传播神经网络分割模型的输入为包含7个特征的特征向量；

所述特征为局部灰度信息和Hu矩不变量；

所述局部信息包括中心像素灰度值、窗口内像素灰度的标准差、中心像素与窗口中最小灰度的绝对差值、中心像素与窗口中最大灰度的绝对差值、以及中心像素与窗口中平均灰度的绝对差值；

所述Hu矩不变量包括|log(I₁)|和|log(I₂)|，其中I₁为每个像素的第一Hu矩，I₂为每个像素的第二Hu矩；

将步骤1.2.3所得血管增强后的图像通过神经网络分割模型进行初步分割的步骤如下：

将步骤1.2.3所得血管增强后的图像输入神经网络分割模型后，神经网络分割模型获取步骤1.2.3所得血管增强后图像的每个像素的特征向量(7个特征)，并以每个像素对应的特征向量(7个特征)作为神经网络分割模型的输入，分别利用神经网络分割模型对眼底图像中每个像素进行二分类(背景像素或视网膜血管像素)，得到神经网络模型的初步分割结果。

作为本发明面向青光眼临床诊断的视网膜血管自动分割方法的进一步改进：

所述步骤2中构建多尺度线检测分割模型的步骤如下：

采用七个尺度s用于直线检测，其算子的长度分别为3、5、7、9、11、13和15；

对步骤1.3预处理后图像中的每个像素，计算以该像素为中心、大小为15×15像素窗口内的平均灰度值I^avg；

在每个尺度s(即上述7个尺度s，取值即为算子的长度)下，对于步骤1.3预处理后图像中的每个像素，在18(如0°、10°、20°...170°)个方向上沿着长度为s的直线(以该像素为中心像素)，计算直线上像素的平均灰度值I；每个像素所有方向的平均灰度值I的最大值作为最大直线检测响应I^s；

每个尺度的直线响应R^s为该尺度的最大直线检测响应I^s与平均灰度值I^avg之间的差值，即R^s＝I^s-I^avg；

对R^s的值进行标准化后得到

，公式如下：

其中，

表示每个尺度所有直线响应的均值，

表示每个尺度所有直线响应的标准差；

注：上述最大响应I^s、平均灰度值I^avg、每个尺度所有直线响应的均值

和标准差

的检测及计算方法为现有技术。

多尺度直线响应R(即多尺度线检测分割模型)：

其中I表示该像素的原始灰度值；

将多尺度直线响应R输出值缩放至0和1之间后进行Otsu法阈值分割，获得多尺度线检测分割模型的初步分割结果。

作为本发明面向青光眼临床诊断的视网膜血管自动分割方法的进一步改进：

所述步骤2中构建尺度空间分析分割模型的步骤如下：

将标准差为δ的高斯核函数G(x,y；δ)与步骤1.1所得绿色分量图像I(x,y)进行卷积，获得每个长度比例因子l的卷积图像I(x,y；l)：

I(x,y；s)＝I(x,y)*G(x,y；δ)

其中，

(σ为高斯核沿x轴坐标中心的离散度)。

对I(x,y；s)求偏导数，得到：

I_x＝I(x,y)*lG_x

I_y＝I(x,y)*lG_y

I_xx＝I(x,y)*l²G_xx

I_xy＝I(x,y)*l²G_xy

I_yy＝I(x,y)*l²G_yy

用梯度的幅值

表示边缘强度，用Hessian矩阵H的最大特征值λ_max计算血管强度，公式如下：

其中，I_xy＝I_yx。

将上述

和H进行归一化处理：

其中l(长度比例因子)的取值范围为1.5≤l≤10，γ表示梯度值，k表示血管强度值；

计算每个像素的血管强度值k，得到视网膜血管的概率图，对其进行Otsu法阈值分割，并均值滤波后得到尺度空间分析分割模型的初步分割结果。

作为本发明面向青光眼临床诊断的视网膜血管自动分割方法的进一步改进：

所述步骤2中构建形态学分割模型的步骤如下：

利用数学形态学中的top-hat变换提取眼底图像中的视网膜血管；

top-hat算子tophat(img)：

tophat(img)＝img-min(open(close(img,s_c),s_o),img)

其中，img为步骤1.1所得绿色分量图像；s_c为形态学“闭”运算的结构元素；open为形态学“开”运算，close为形态学“闭”运算；

所述形态学“闭”运算(即，close)用于消除因小噪声波动引起的影响；

所述形态学“开”运算(即，open)的结构元素分别取半径为1、2、3、4、5、6、7、8像素的圆盘，获得8幅图像，之后对结构元素为连续半径的两幅图像求平均值，将其缩减到4幅图像；

之后令每个像素的响应除以对应尺度值的平方，所有尺度的最大响应值作为该像素的最终响应，并将其归一化至[0,1]区间获得形态学分割模型的初步分割结果。

作为本发明面向青光眼临床诊断的视网膜血管自动分割方法的进一步改进：

所述步骤3中创建掩膜的方法为：

按照步骤1.1提取待分割眼底图像的绿色分量图像，并对待分割眼底图像的绿色分量图像进行自适应直方图均衡化；

对自适应直方图均衡化后的图像进行阈值分割得到掩膜图像(阈值τ_mask设置为0.73)，并使用圆盘状结构元素(半径为3像素)对掩膜图像的白色区域进行形态学“膨胀”操作。

作为本发明面向青光眼临床诊断的视网膜血管自动分割方法的进一步改进：

所述步骤3中多模型融合的方法为：

形态学分割模型的初步分割结果和经过“膨胀”操作后的掩膜图像进行逻辑“与”运算，掩膜图像中非白色区域和步骤2中初始分割输出进行逻辑“与”运算，两个运算结果进行逻辑“或”运算后得到多模型融合的视网膜血管分割的组合结果；

作为本发明面向青光眼临床诊断的视网膜血管自动分割方法的进一步改进：

所述步骤4中精分割方法包括以下步骤：

4.1、种子像素的选取：利用Otsu法将组合结果中的像素分为视网膜血管和背景两类；所有灰度值小于阈值t的像素归为初始背景类；反之，所有强度值大于或等于阈值t的像素则归为初始视网膜血管类；

利用灰度均值μ_b和标准差σ_b表征背景类像素，灰度均值μ_v和标准差σ_v表征背景类表征视网膜血管类像素；若像素的灰度值大于μ_v，则认为该像素是视网膜血管类像素；反之，若像素的强度值小于μ_v，则认为该像素是背景类像素；将这些具有确定类别标签的像素作为区域生长法的种子像素点；

4.2、迭代生长：利用区域生长法和梯度信息对背景和视网膜血管区域进行迭代生长和分类；

对于血管类像素：

对于背景类像素：

其中，

为像素点的灰度值；γ为像素点的梯度值；

初始值α＝1，每迭代生长一次，其值增加0.6，直至所有像素分类完毕；

然后于不考虑梯度信息的条件下再次进行区域迭代生长，直至所有像素均分类完毕，此时得到视网膜血管的精分割结果。

本发明与现有技术相比，具有如下技术优势：

1、本发明融合了匹配滤波器模型、神经网络模型、多尺度线检测模型、尺度空间分析模型和形态学模型五种依赖于不同图像处理技术的模型，很好地消除了视盘和渗出物等明亮区域的影响，能获得理想的分割结果；

2、本发明不需要海量级别的数据建立视网膜血管分割模型，极大地降低了需要处理的数据量和复杂度，易于实现，能有效提高视网膜血管分割的效率，满足临床实际应用对实时性的要求；

3、本发明考虑视网膜血管的属性(即血管网络是由连通血管段的血管树构成)，在多模态融合结果的基础上利用区域生长法和梯度信息对背景和血管区域进行迭代生长，分割结果具有较好的连续性和平滑性，能保留更多的视网膜血管细节和更完整视的网膜血管网络。

附图说明

下面结合附图对本发明的具体实施方式作进一步详细说明。

图1为本发明面向青光眼临床诊断的视网膜血管自动分割方法的原理流程图；

图2为本发明实施例中神经网络分割模型的预处理结果(图中a为绿色分量图像，b为背景图像，c为背景均匀化后的图像，d为血管增强后的图像)；

图3为本发明实施例中匹配滤波器分割模型的视网膜血管分割结果；

图4为本发明实施例中神经网络模型的视网膜血管分割结果；

图5为本发明实施例中多尺度线检测模型的视网膜血管分割结果；

图6为本发明实施例中尺度空间分析分割模型的视网膜血管分割结果；

图7为本发明实施例中形态学模型的视网膜血管分割结果；

图8为本发明实施例中多模型融合的视网膜血管组合结果(图中a形态学模型分割结果和掩膜白色区域进行逻辑“与”运算后的图像，b为初步分割输出，c为初步分割输出和掩膜中其他区域进行逻辑“与”运算后的图像，d为多模型融合的组合结果)；

图9为本发明实施例中视网膜血管的精分割结果；

图10为本发明与不同方法对DRIVE眼底图像库测试集中编号为03_test、06_test和18_test三个示例图像进行视网膜血管分割的可视化结果对比。

具体实施方式

下面结合具体实施例对本发明进行进一步描述，但本发明的保护范围并不仅限于此。

实施例1、面向青光眼临床诊断的视网膜血管自动分割方法，如图1-9所示，包括以下内容：

本发明首先将眼底图像预处理，之后分别构建匹配滤波器、神经网络、多尺度线检测、尺度空间分析和形态学模型初步分割视网膜血管，为减少噪声取五个分割结果的均值作为初步分割输出。其次设计掩膜分离渗出物和视盘区域，将形态学模型分割结果替换掩膜的白色区域，并融合初步分割输出生成组合结果。最后考虑视网膜血管先验知识(即，视网膜血管网络是由连通血管段的血管树构成)，对组合结果阈值分割和区域迭代生长后获取最终结果。

本实施例采用的测试数据为国际上公认的DRIVE眼底图像库测试集(均为RGB格式的24位真彩色眼底图像)中编号为02_test的眼底图像。

注：DRIVE眼底图像库总共含有40幅彩色眼底图像，其中7幅是早期轻度DR眼底图像。该眼底图像库的所有图像的分辨率为768×584像素，分成训练图像集和测试图像集，各含有20幅眼底图像。对于训练图像集，DRIVE提供了一个单一的手动分割结果作为真值图；对于测试图像集，则提供了manual1和manual2(2nd observer)两种手动分割结果。

注：实际待分割眼底图像的处理方式等同于本实施例中对该测试数据的处理方式，故本法不再对实际待分割眼底图像的处理方式进行详细描述。

具体工作内容包括如下步骤：

步骤1、眼底图像预处理：分别设计适合匹配滤波器模型、神经网络模型、多尺度线检测模型、尺度空间分析模型和形态学模型的眼底图像预处理步骤，以降低噪声影响。本发明针对匹配滤波器模型、神经网络模型、多尺度线检测模型、尺度空间分析模型和形态学模型采取不同的预处理步骤，具体包括以下步骤：

1.1、匹配滤波器模型、尺度空间分析模型以及形态学分割模型的预处理：

将彩色的眼底图像分解为红色、绿色和蓝色三个分量图像，并提取该眼底图像的绿色分量图像。

注：由于绿色分量图像中视网膜血管与背景的对比度最高，因此提取眼底图像的绿色分量图像进行视网膜血管分割。

1.2、神经网络分割模型的预处理：

按照步骤1.1中提取眼底图像的绿色分量图像，为了减轻中央血管光反射的影响，还对该绿色分量图像按照以下3步进行预处理：

1.2.1、利用直径为3像素的圆盘状结构元素对绿色分量图像进行形态学“开”运算，并对所得图像进行均值滤波，窗口大小为69×69像素，得到背景图像；

1.2.2、利用按照步骤1.1所提取的绿色分量图像减去步骤1.2.1所得的背景图像，并将像素灰度值调整至0-1范围内实现背景均匀化，在整个图像集上生成统一的背景灰度级，获得背景均匀化后的图像。

1.2.3、首先将步骤1.2.2所得背景均匀化后的图像中所有像素的灰度值加上常数0.5，使图像中血管增强，之后对该图像进行“补”运算，并使用半径为8像素的圆盘状结构元素对“补”运算后的图像进行top-hat(顶帽)变换，从而增强该图像的暗区域以及视网膜血管区域，同时去除明亮区域(如视盘区域)，获得血管增强后的图像。

预处理结果如图2所示，图2a为绿色分量图像；图2b为背景图像；图2c为背景均匀化后的图像；图2d为血管增强后的图像。

1.3、多尺度线检测分割模型的预处理：

首先按照步骤1.1中提取眼底图像的绿色分量图像，并对该绿色分量图像进行“取反”操作，之后对所得图像进行背景均匀化处理。之后将背景均匀化处理后的图像进行均值滤波(窗口大小为69×69像素)，并将均值滤波后的局部平均灰度替换明亮区域，以限制视盘的影响，获得预处理后的眼底图像。

注：实现背景均匀化的方法如步骤1.2.1所示。

步骤2：初步分割：分别构建匹配滤波器模型、神经网络模型、多尺度线检测模型、尺度空间分析模型和形态学模型，初步分割视网膜血管，得到各模型的初步分割结果；为了减少噪声，取五个模型分割结果的均值作为初步分割输出。

具体包括以下步骤：

2.1、构建匹配滤波器分割模型：

本发明所采用匹配滤波器分割模型为二维匹配滤波器。

在眼底图像中，视网膜血管正面的横截面的亮度分布可近似为高斯函数，而高斯曲线的表达式为：

其中，K(x,y)称为高斯核函数(下文中简称为高斯核)；(x,y)表示图像(步骤1.1所得绿色分量图像)像素的坐标，满足|x|≤λσ，|y|≤L/2；σ为高斯核沿x轴坐标中心的离散度，取值为2；L为沿y轴被截断视网膜血管的长度，取值为9；λ为常数，本实施例设置为3。由于视网膜血管的方向是任意的，故高斯核应旋转以匹配不同方向的视网膜血管。

二维匹配滤波器的高斯核是一个11×11的旋转矩阵(即，高斯核旋转矩阵)，假设g＝[x,y]是高斯核旋转矩阵中的离散点，θ_i表示第i个高斯核的旋转角度。在计算高斯核旋转矩阵中的系数时，假设旋转中心在(0,0),则旋转矩阵r_i为：

每隔15°就选取一个高斯核，这样在0°～180°范围内总共选取了12个不同的高斯核旋转矩阵。

将步骤1.1所得绿色分量图像通过匹配滤波器分割模型进行初步分割的步骤如下：

将步骤1.1预处理后的眼底图像(即，绿色分量图像)中每个像素分别与12个不同的高斯核旋转矩阵进行卷积，选取其中最大的卷积值作为匹配滤波后的像素值，得到视网膜血管概率图；之后对视网膜血管概率图进行阈值分割(全局阈值t_mf设置为0.33)，得到含有不同区域(即，背景区域与非背景区域)的二值化图像。若该像素所处区域面积大于950像素，且该区域为非背景区域，则将该像素标记为视网膜血管像素进行合并，即，将所处区域面积大于950像素，且该区域为非背景区域的像素标记为视网膜血管像素，合并所有视网膜血管像素，获得匹配滤波器分割模型的初步分割结果。

图3为σ＝2，L＝9，λ＝3时，匹配滤波器分割模型对本实施例中步骤1.1提取的眼底图像的绿色分量图像进行视网膜血管分割所得结果。

匹配滤波器分割模型很好地利用了视网膜血管横截面与高斯曲线吻合的特性，在一定程度上能够平滑噪声，使得视网膜血管结构得到增强。

2.2、构建神经网络分割模型：

本发明神经网络分割模型为多层的前馈反向传播神经网络分割模型，包含1个输入层(7个节点)、3个隐层(15个节点)和1个输出层(1个节点)。隐层的传输函数采用线性函数，输出层的传输函数logsig(z)采用对数S型函数：

其中，z表示隐层的输出值。

神经网络分割模型的输入为包含7个特征的特征向量，分别为5个局部灰度信息和2个Hu矩不变量。5个局部灰度信息包括中心像素灰度值、窗口内像素灰度的标准差、中心像素与窗口中最小灰度的绝对差值、中心像素与窗口中最大灰度的绝对差值、以及中心像素与窗口中平均灰度的绝对差值。

注：对于位置为(x,y)处的像素，其灰度特征通过计算以该像素为中心、大小为9×9像素窗口内的灰度获得。

上述5个局部灰度信息均可通过灰度特征均通过现有技术计算获得。

对于每个像素，其第一Hu矩和第二Hu矩计算步骤如下：

若血管增强后的图像区域为f(x,y)，其(p+q)阶原点矩m_pq和(p+q)中心矩μ_pq分别为:

其中，p,q＝0,1,2,…；

为目标区域灰度质心，

f(x,y)归一化后的(p+q)中心矩η_pq为：

其中，

信息量最大的前两个Hu矩I₁和I₂分别为：

I₁＝η₂₀+η₀₂

本实施例为了使2个Hu矩不变量具有零均值和单位方差，取Hu矩I₁和I₂对数的绝对值|log(I₁)|和|log(I₂)|作为位置为(x,y)处像素最终的两个特征(即，Hu矩不变量)。

将步骤1.2.3所得血管增强后的图像通过神经网络分割模型进行初步分割的步骤如下：

将步骤1.2.3所得血管增强后的图像输入神经网络分割模型后，神经网络分割模型根据上述步骤，获取步骤1.2.3所得血管增强后图像的每个像素的特征向量(7个特征)，并以每个像素对应的特征向量(7个特征)作为神经网络分割模型的输入，分别利用神经网络分割模型对眼底图像中每个像素进行二分类(背景像素或视网膜血管像素)，从而实现对视网膜血管的输出。

利用构建的神经网络分割模型对步骤1.2.3所得血管增强后的图像进行视网膜血管分割，所得结果如图4所示。神经网络分割模型具有很强的非线性拟合能力和鲁棒性，提取的视网膜血管网络连通性较好。

2.3、构建多尺度线检测分割模型：

采用七个尺度s用于直线检测，其算子的长度分别为3、5、7、9、11、13和15像素。具体处理步骤如下：

对于步骤1.3预处理后图像中的每个像素，计算以该像素为中心、大小为15×15像素窗口内的平均灰度值I^avg。在每个尺度s(即上述7个尺度s，取值即为算子的长度)下，对于步骤1.3预处理后图像中的每个像素，在18(0°、10°、20°...170°)个方向上沿着长度为s的直线(以该像素为中心像素)，计算直线上像素的平均灰度值

每个像素所有方向的平均灰度值

的最大值作为最大直线检测响应I^s。每个尺度的直线响应R^s为该尺度的最大直线检测响应I^s与平均灰度值I^avg之间的差值，即R^s＝I^s-I^avg。为了使R^s具有零均值和单位方差，本发明对R^s的值进行标准化后得到

其中，

表示每个尺度所有直线响应的均值，

表示每个尺度所有直线响应的标准差。

注：上述最大响应I^s(x,y)、平均灰度值I^avg、每个尺度所有直线响应的均值

和标准差

的检测及计算方法为现有技术，故不在本说明书中进行详细介绍。

在完成直线响应

的计算后，多尺度直线响应R(即多尺度线检测分割模型)就可以通过每个尺度s的直线响应

和图像中原始灰度值I(即像素的原始灰度值)的线性组合获得：

将步骤1.3预处理后的眼底图像通过多尺度线检测分割模型进行初步分割的步骤如下：

对步骤1.3预处理后的眼底图像的每个像素，计算以该像素为中心、大小为15×15像素窗口内的平均灰度值I^avg；对于每个尺度s，在18个方向上沿着以该像素为中心像素、长度为s的直线，计算直线上像素的平均灰度值

其最大值作为最大直线检测响应I^s；将该尺度的最大直线检测响应I^s减去平均灰度值I^avg获得每个尺度的直线响应R^s，并对R^s的值进行标准化后得到

通过每个尺度s的直线响应

和图像中原始灰度值I(即像素的原始灰度值)的线性组合获得多尺度线线检测分割模型R。将多尺度线线检测分割模型最终输出值缩放至0和1之间后进行Otsu法阈值分割(即，灰度图像的自动阈值分割，为现有技术，故不详细介绍)，以保留图像中更精细的特征，并移除面积小于100像素的所有连通区域。

利用构建的多尺度线检测分割模型对步骤1.3预处理后的眼底图像进行视网膜血管分割，所得结果如图5所示。多尺度线检测分割模型能够消除视网膜血管中间区域的光反射，避免紧邻视网膜血管之间互相干扰。

2.4、构建尺度空间分析分割模型：

对于步骤1.1预处理所得绿色分量图像I(x,y)，将标准差为δ的高斯核函数G(x,y；δ)与步骤1.1所得绿色分量图像I(x,y)进行卷积，获得每个长度比例因子l的卷积图像I(x,y；l)：

I(x,y；l)＝I(x,y)*G(x,y；δ)

其中，

为了计算梯度

和Hessian矩阵H，对I(x,y；l)求偏导数，得到：

I_x＝I(x,y)*lG_x

I_y＝I(x,y)*lG_y

I_xx＝I(x,y)*l²G_xx

I_xy＝I(x,y)*l²G_xy

I_yy＝I(x,y)*l²G_yy

若边缘强度用梯度的幅值

表示，那么血管的强度可通过Hessian矩阵H的最大特征值λ_max估算。即：

其中，I_xy＝I_yx。

注：上述计算Hessian矩阵H的最大特征值λ_max为现有技术，故不在本说明书中进行详细介绍。

由于大血管的局部最大响应比小血管的局部最大响应要大得多，而长度比例因子l与视网膜血管的尺寸有关，故本发明将边缘强度和血管强度进行归一化，即：

完成上述步骤之后，本发明根据图像中最小血管和最大血管的尺寸设定长度比例因子l的取值范围(1.5≤l≤10)，以0.5为步长在区间[1.5,10]内计算尺度空间信息γ和k，并获取γ和k的最大值；γ最大值作为步骤4.2中区域迭代生长的梯度信息，k最大值作为该像素点的最终强度值。

将步骤1.1预处理后的眼底图像通过尺度空间分析分割模型进行初步分割的步骤如下：

利用高斯核函数对步骤1.1预处理所得绿色分量图像进行卷积，获取每个长度比例因子s的卷积图像；计算Hessian矩阵，得到其最大特征值作为血管强度的估算，并进行归一化处理；计算尺度空间信息，并获取其最大值作为该像素点的最终强度值。对每个像素进行相同操作后，得到视网膜血管的概率图，对其进行Otsu法阈值分割，并进行窗口大小为3×3像素的均值滤波，得到视网膜血管分割结果。

利用构建的尺度空间分析分割模型对步骤1.1预处理后的眼底图像进行视网膜血管分割，所得结果如图6所示。尺度空间分析分割模型能够实现视网膜血管增强，并适应视网膜血管宽度的变化。

2.5、构建形态学分割模型：

利用数学形态学中的top-hat变换提取步骤1.1所得绿色分量图像中的视网膜血管。由于常规的top-hat算子对噪声异常敏感，故采用修改后的top-hat算子tophat(img)构建形态学分割模型：

tophat(img)＝img-min(open(close(img,s_c),s_o),img)

其中，img为步骤1.1所得绿色分量图像；s_c为形态学“闭”运算的结构元素；s_o为形态学“开”运算的结构元素。

形态学“闭”运算(即，close)能够消除因小噪声波动引起的影响，其结构元素尺寸设定为比小噪声稍大，本实施例中小噪音尺寸为1像素，其结构元素尺寸为半径为2像素的圆盘；

形态学“开”运算(即，open)结构元素分别取半径为1、2、3、4、5、6、7、8像素的圆盘，从而获得8幅不同尺度的图像。

本发明对结构元素为连续半径的两幅图像求平均值，这样不同尺度的8幅图像最终缩减到4幅图像(即尺度为1、2、3和4)，有助于减少噪声。为了弥补较弱的响应，每个像素的响应除以对应尺度值的平方，所有尺度的最大响应值作为该像素的最终响应，并将其归一化至[0,1]区间即得到输出结果。

将步骤步骤1.1所得绿色分量图像通过形态学分割模型进行初步分割的步骤如下：

首先对步骤1.1所得绿色分量图像进行形态学“闭”运算(即，close)，消除因小噪声波动引起的影响；其次将完成形态学“开”运算的图像进行形态学“开”运算(即，open)，其结构元素分别取半径为1、2、3、4、5、6、7、8像素的圆盘，从而获得8幅图像。然后为了进一步减少噪声，对结构元素为连续半径的两幅图像求平均值，得到4幅不同尺度的图像(尺度为1、2、3和4)。为了弥补较弱的响应，每个像素的响应除以对应尺度值的平方，所有尺度的最大响应值作为该像素的最终响应，并将其归一化至[0,1]区间即得到输出结果。

本实施例利用构建的形态学分割模型对步骤1.1预处理后的眼底图像进行视网膜血管分割，所得结果如图7所示。形态学分割模型运算速度快，抗噪能力强。

2.6、计算匹配滤波器模型、神经网络模型、多尺度线检测模型、尺度空间分析模型和形态学模型五种分割模型的视网膜血管初步分割结果的均值，得到初步分割输出。

3、多模型融合：为发挥上述五个模型的各自优点，本发明还设计掩膜，分离眼底图像中渗出物和视盘区域，并利用形态学模型分割结果替换掩膜中白色区域，最后将所得图像与初步分割输出融合，生成组合结果。具体步骤如下：

3.1、创建掩膜：

按照步骤1.1提取眼底图像的绿色分量图像；为了消除因不均匀光照引起像素灰度变化的影响，对眼底图像的绿色分量图像进行自适应直方图均衡化(直方图均衡化为现有技术，故不详细介绍)；接着对自适应直方图均衡化后的图像进行阈值分割得到掩膜图像，本实施例中阈值τ_mask设置为0.73；同时，为了消除视盘区域边缘效应的影响，使用半径为3像素的圆盘状结构元素对掩膜图像的白色区域进行形态学“膨胀”操作。

3.2、多模型融合：

将步骤2.5中形态学模型的视网膜血管分割结果和步骤3.1中白色区域“膨胀”操作后的掩膜图像进行逻辑“与”运算，掩膜图像中其他区域(即，非白色区域)和步骤2.6中初始分割输出进行逻辑“与”运算，两个运算结果进行逻辑“或”运算后得到多模型融合的视网膜血管分割的组合结果。如图8所示，图8a为形态学模型分割结果和掩膜白色区域进行逻辑“与”运算后的图像，图8b为初步分割输出，图8c为初步分割输出和掩膜中其他区域进行逻辑“与”运算后的图像，图8d为多模型融合的组合结果。

4、精分割：考虑视网膜血管先验知识，应用Otsu法阈值分割组合结果，并根据血管的连通特性进行区域迭代生长后获取视网膜血管分割的最终结果。包括以下具体步骤：

4.1、种子像素的选取：

利用Otsu法将组合结果中的像素分为视网膜血管和背景两类。所有灰度值小于阈值t的像素归为初始背景类；反之，所有强度值大于或等于阈值t的像素则归为初始视网膜血管类。利用各自的灰度均值μ_b和μ_v以及标准差σ_b和σ_v分别表征背景类和视网膜血管类像素。若像素的灰度值大于μ_v，则认为该像素是视网膜血管类像素；反之，若像素的强度值小于μ_v，则认为该像素是背景类像素。将这些具有确定类别标签的像素作为区域生长法的种子像素点。

注：灰度均值μ_b和标准差σ_b表征背景类像素，灰度均值μ_v和标准差σ_v表征背景类表征视网膜血管类像素。灰度均值μ_b和μ_v以及标准差σ_b和σ_v均通过现有技术计算获得。

4.2、迭代生长：

利用区域生长法和梯度信息对背景和视网膜血管区域进行迭代生长和分类。

由于梯度低的像素通常在视网膜血管或者背景的中间区域，因此这些像素更易分类。本文利用均值μ_g和标准差σ_g表征梯度值的直方图。具体的迭代生长过程如下：

首先，本文只对那些具有明确概率值(非常高或非常低)和低梯度值的像素进行分类。

即：

对于血管类像素：

对于背景类像素：

其中，

为像素点的灰度值；γ为像素点的梯度值。本文取初始值α＝1，每生长一次，其值增加0.6，直至所有像素分类完毕。视网膜血管外的边缘像素通常具有高梯度值，而条件γ≤μ_g则可防止将这些边缘像素归为背景类。均值μ_g和标准差σ_g通过现有技术计算获得。

然后，再次进行区域迭代生长，但不考虑梯度信息，直至所有像素均分类完毕，此时得到视网膜血管的精分割结果。图9为本发明取t＝0.1667、μ_b＝0.0087、μ_v＝0.3258、μ_g＝0.0789、σ_b＝0.0211、σ_v＝0.0908和σ_g＝0.0774时，在步骤3.2中组合结果的基础上进行区域迭代生长得到的视网膜血管精分割结果。

实际待分割眼底图像的分割步骤与上述对编号为02_test的眼底图像的分割步骤完全相同，故不在本说明书中进行重复告知。

实验1：定量分析和检验不同视网膜血管分割算法的分割效果；

将实施例1所得分割结果与DRIVE眼底图像库所提供的manual1手动分割结果进行比较，分析像素的分类情况。本发明采用常用的检测精度A_cc(Accuracy)、敏感度S_n(Sensitivity)和特异性S_p(Specificity)三个指标客观评价视网膜血管分割算法的性能，分别为：

其中，TP(true positive)为真阳性，即视网膜血管像素点分割正确的个数；FP(false positive)为假阳性，即非视网膜血管像素点被错误分割为视网膜血管像素点的个数；TN(true negative)为真阴性，即非视网膜血管像素点分割正确的个数；FN(falsenegative)为假阴性，即视网膜血管像素点被错误分割为非视网膜血管像素点的个数。结果如表1所示：

按照现有非监督算法2nd observer、Chaudhuri、Mendonca、Zhang、Fraz和Zhao分别对DRIVE眼底图像库测试集进行视网膜血管分割所得检测精度A_cc、敏感度S_n和特异性S_p结果的定量对比，结果如表1所示。

上述Zhang为基于匹配滤波的视网膜血管分割算法，Fraz为基于形态学的视网膜血管分割算法，Zhao为基于形变模型的视网膜血管分割算法

表1

由表1可以看出，本发明分割视网膜血管的准确性较高，与其他经典算法相比具有一定的竞争力，可望为临床决策提供可靠的参考。其中，在检测精度A_cc方面，本发明性能要优于Chaudhuri、Zhang、Fraz三种非监督算法和Staal监督算法，但要略逊于Mendonca、Zhao两种非监督算法和Zhu、Soares等监督算法；在敏感度S_n方面，本发明则表现出众，优于Chaudhuri、Mendonca、Zhao等非监督算法和Zhu、Franklin、Soares等监督算法，这说明本文算法对视网膜血管的生长有效。此外，本发明的特异性S_p要比Chaudhuri、Mendonca、Zhang等非监算法以及Soares、Wang等监督算法高，略低于Zhu和Franklin监督算法。

实验2：分割后视网膜血管网络可视化对比情况：

分别利用本发明即现有分割方法对DRIVE眼底图像库测试集中编号为03_test、06_test和18_test三个示例图像进行视网膜血管分割，并将所得结果进行可视化对比。

注：上述现有分割方法为：Chaudhuri算法、Jiang算法、Niemeijer算法和staal算法。

由图10可知，本发明分割眼底图像所得视网膜血管网络的完整度较高，不仅视网膜血管的主干和末梢的连通性好，而且较好地分割出大多数微细视网膜血管。

本发明融合了匹配滤波器、神经网络、多尺度线检测、尺度空间分析和形态学五种依赖于不同图像处理技术的模型，很好地消除了视盘和渗出物等明亮区域的影响，能获得理想的分割结果。同时，本发明不需要海量级别的数据建立视网膜血管分割模型，极大地降低了需要处理的数据量和复杂度，易于实现，能有效提高视网膜血管分割的效率。此外，考虑视网膜血管的属性(即血管网络是由连通血管段的血管树构成)，本发明在多模态融合结果的基础上利用区域生长法和梯度信息对背景和血管区域进行迭代生长，分割结果具有较好的连续性和平滑性，能保留更多的视网膜血管细节和更完整视的网膜血管网络。视网膜血管的自动分割技术不但能减轻眼科医生的负担，而且能有效解决偏远地区缺乏有经验眼科医生的问题。因此，视网膜血管自动分割技术对于青光眼的辅助诊疗具有极其重要的意义。

最后，还需要注意的是，以上列举的仅是本发明的若干个具体实施例。显然，本发明不限于以上实施例，还可以有许多变形。本领域的普通技术人员能从本发明公开的内容直接导出或联想到的所有变形，均应认为是本发明的保护范围。

Claims (8)

1.面向青光眼的视网膜血管自动分割方法，其特征在于包括以下步骤：

步骤1：眼底图像预处理：预处理待分割眼底图像；

眼底图像预处理针对匹配滤波器模型、神经网络模型、多尺度线检测模型、尺度空间分析模型和形态学模型采取不同预处理步骤：

1.1、匹配滤波器模型、尺度空间分析模型以及形态学分割模型的预处理：

将待分割眼底图像分解为红色、绿色和蓝色三个分量图像，提取待分割眼底图像的绿色分量图像；

1.2、神经网络分割模型的预处理：

按照步骤1.1提取待分割眼底图像的绿色分量图像，并对其进行以下处理：

1.2.1、利用圆盘状结构元素对绿色分量图像进行形态学“开”运算，并对所得图像进行均值滤波，得到背景图像；

1.2.2、利用按照步骤1.1所提取的绿色分量图像减去步骤1.2.1所得的背景图像，并将像素灰度值调整至0-1范围内实现背景均匀化，获得背景均匀化后的图像；

1.2.3、首先将步骤1.2.2所得背景均匀化后的图像中血管进行增强，之后对该图像进行“补”运算，并使用圆盘状结构元素对“补”运算后的图像进行top-hat变换，增强该图像的暗区域以及视网膜血管区域，同时去除明亮区域，获得血管增强后的图像；

1.3、多尺度线检测分割模型的预处理：

首先按照步骤1.1中提取眼底图像的绿色分量图像，并对该绿色分量图像进行“取反”操作，之后对所得图像进行背景均匀化处理；之后将处理后的图像进行均值滤波，并将均值滤波后的局部平均灰度替换明亮区域；

步骤2：初步分割：分别构建匹配滤波器模型、神经网络模型、多尺度线检测模型、尺度空间分析模型和形态学模型；利用上述构建的模型，分别对步骤1所得待分割眼底图像进行初步分割，获得对应的初步分割结果；将五个初步分割结果的均值作为初步分割输出；

步骤3：多模型融合：利用掩膜分离步骤1形态学模型对应预处理后的待分割眼底图像中渗出物和视盘区域，并利用步骤2所得的形态学模型分割结果替换掩膜中白色区域，将掩膜图像中的非白色区域和初始分割输出进行逻辑“与”运算，再将两个运算结果进行逻辑“或”运算后得到多模型融合的视网膜血管分割的组合结果；

步骤4：精分割：利用Otsu法阈值分割步骤3所得组合结果，并根据血管的连通特性进行区域迭代生长后，获取视网膜血管分割的最终结果。

2.根据权利要求1所述的面向青光眼的视网膜血管自动分割方法，其特征在于：

所述步骤2中构建匹配滤波器分割模型步骤如下：

所述匹配滤波器分割模型为二维匹配滤波器，其高斯曲线的表达式为：

其中，K(x,y)称为高斯核函数；(x,y)表示步骤1.1所得绿色分量图像中像素的坐标，满足|x|≤λσ，|y|≤L/2；σ为高斯核沿x轴坐标中心的离散度；L为沿y轴被截断视网膜血管的长度；λ为常数；

所述高斯核函数为旋转矩阵，在0°～180°范围内，每隔15°就选取一个高斯核函数，共选取了12个不同的高斯核旋转矩阵；

将步骤1.1所得绿色分量图像通过匹配滤波器分割模型进行分割的步骤如下：

将步骤1.1高斯核旋转矩阵中每个像素分别与12个不同的高斯核旋转矩阵进行卷积，选取其中最大的卷积值作为匹配滤波后的像素值，得到视网膜血管概率图；之后对视网膜血管概率图进行阈值分割，得到含有背景区域及非背景区域的二值化图像；

将所处区域面积大于950像素，且该区域为非背景区域的像素标记为视网膜血管像素，合并所有视网膜血管像素，获得匹配滤波器分割模型的初步分割结果。

3.根据权利要求1所述的面向青光眼的视网膜血管自动分割方法，其特征在于：

所述步骤2中构建神经网络分割模型的步骤如下：

前馈反向传播神经网络为多层的前馈反向传播神经网络分割模型，包含1个输入层、3个隐层和1个输出层；隐层的传输函数采用线性函数，输出层的传输函数logsig(z)采用对数S型函数：

其中，z表示隐层的输出值；

所述前馈反向传播神经网络分割模型的输入为包含7个特征的特征向量；

所述特征为局部灰度信息和Hu矩不变量；

局部信息包括中心像素灰度值、窗口内像素灰度的标准差、中心像素与窗口中最小灰度的绝对差值、中心像素与窗口中最大灰度的绝对差值、以及中心像素与窗口中平均灰度的绝对差值；

所述Hu矩不变量包括|log(I₁)|和|log(I₂)|，其中I₁为每个像素的第一Hu矩，I₂为每个像素的第二Hu矩；

将步骤1.2.3所得血管增强后的图像通过神经网络分割模型进行初步分割的步骤如下：

将步骤1.2.3所得血管增强后的图像输入神经网络分割模型后，神经网络分割模型获取步骤1.2.3所得血管增强后图像的每个像素的特征向量，并以每个像素对应的特征向量作为神经网络分割模型的输入，分别利用神经网络分割模型对眼底图像中每个像素进行二分类，得到神经网络模型的初步分割结果。

4.根据权利要求1所述的面向青光眼的视网膜血管自动分割方法，其特征在于：

所述步骤2中构建多尺度线检测分割模型的步骤如下：

采用七个尺度s用于直线检测，其算子的长度分别为3、5、7、9、11、13和15；

对步骤1.3预处理后图像中的每个像素，计算以该像素为中心、大小为15×15像素窗口内的平均灰度值I^avg；

在每个尺度s下，对于步骤1.3预处理后图像中的每个像素，在18个方向上沿着长度为s的直线，计算直线上像素的平均灰度值

每个像素所有方向的平均灰度值I的最大值作为最大直线检测响应I^s；

每个尺度的直线响应R^s为该尺度的最大直线检测响应I^s与平均灰度值I^avg之间的差值，即R^s＝I^s-I^avg；

对R^s的值进行标准化后得到

公式如下：

其中，

表示每个尺度所有直线响应的均值，

表示每个尺度所有直线响应的标准差；

多尺度直线响应R：

其中，I表示该像素的原始灰度值；

将多尺度直线响应R输出值缩放至0和1之间后进行Otsu法阈值分割，获得多尺度线检测分割模型的初步分割结果。

5.根据权利要求1所述的面向青光眼的视网膜血管自动分割方法，其特征在于：

所述步骤2中构建尺度空间分析分割模型的步骤如下：

将标准差为δ的高斯核函数G(x,y；δ)与步骤1.1所得绿色分量图像I(x,y)进行卷积，获得每个长度比例因子l的卷积图像I(x,y；l)：

I(x,y；l)＝I(x,y)*G(x,y；δ)

其中，

对I(x,y；l)求偏导数，得到：

I_x＝I(x,y)*lG_x

I_y＝I(x,y)*lG_y

I_xx＝I(x,y)*l²G_xx

I_xy＝I(x,y)*l²G_xy

I_yy＝I(x,y)*l²G_yy

用梯度的幅值

表示边缘强度，用Hessian矩阵H的最大特征值λ_max计算血管强度，公式如下：

其中，I_xy＝I_yx；

将上述

和H进行归一化处理：

其中l的取值范围为1.5≤l≤10，γ表示梯度值，k表示血管强度值；

计算每个像素的血管强度值k，得到视网膜血管的概率图，对其进行Otsu法阈值分割，并均值滤波后得到尺度空间分析分割模型的初步分割结果。

6.根据权利要求1所述的面向青光眼的视网膜血管自动分割方法，其特征在于：

所述步骤2中构建形态学分割模型的步骤如下：

利用数学形态学中的top-hat变换提取眼底图像中的视网膜血管；

top-hat算子tophat(img)：

tophat(img)＝img-min(open(close(img,s_c),s_o),img)

其中，img为步骤1.1所得绿色分量图像；s_c为形态学“闭”运算的结构元素；open为形态学“开”运算，close为形态学“闭”运算；

所述形态学“闭”运算用于消除因小噪声波动引起的影响；

所述形态学“开”运算的结构元素分别取半径为1、2、3、4、5、6、7、8像素的圆盘，获得8幅图像，之后对结构元素为连续半径的两幅图像求平均值，将其缩减到4幅图像；

之后令每个像素的响应除以对应尺度值的平方，所有尺度的最大响应值作为该像素的最终响应，并将其归一化至[0,1]区间获得形态学分割模型的初步分割结果。

7.根据权利要求1～6任一所述的面向青光眼的视网膜血管自动分割方法，其特征在于：

所述步骤3中创建掩膜的方法为：

按照步骤1.1提取待分割眼底图像的绿色分量图像，并对待分割眼底图像的绿色分量图像进行自适应直方图均衡化；

对自适应直方图均衡化后的图像进行阈值分割得到掩膜图，并使用圆盘状结构元素对掩膜图像的白色区域进行形态学“膨胀”操作。

8.根据权利要求7所述的面向青光眼的视网膜血管自动分割方法，其特征在于：

所述步骤4中精分割方法包括以下步骤：

4.1、种子像素的选取：利用Otsu法将组合结果中的像素分为视网膜血管和背景两类；所有灰度值小于阈值t的像素归为初始背景类；反之，所有强度值大于或等于阈值t的像素则归为初始视网膜血管类；

利用灰度均值μ_b和标准差σ_b表征背景类像素，灰度均值μ_v和标准差σ_v表征背景类表征视网膜血管类像素；若像素的灰度值大于μ_v，则认为该像素是视网膜血管类像素；反之，若像素的强度值小于μ_v，则认为该像素是背景类像素；将这些具有确定类别标签的像素作为区域生长法的种子像素点；

4.2、迭代生长：利用区域生长法和梯度信息对背景和视网膜血管区域进行迭代生长和分类；

对于血管类像素：

对于背景类像素：

其中，

为像素点的灰度值；γ为像素点的梯度值；

初始值α＝1，每迭代生长一次，其值增加0.6，直至所有像素分类完毕；

然后于不考虑梯度信息的条件下再次进行区域迭代生长，直至所有像素均分类完毕，此时得到视网膜血管的精分割结果。

Images