CN112465842B - 基于U-net网络的多通道视网膜血管图像分割方法 - Google Patents

Abstract

本发明公开了一种基于U‑net网络的多通道视网膜血管图像分割方法。本发明首先对数据集图像进行扩增处理和一系列的预处理，提高图像质量；其次将多尺度匹配滤波算法与改进的形态学算法相结合，构建U‑net网络的多通道特征提取结构；随后对三个通道进行网络训练，获得所需的分割网络，并对输出结果采用自适应阈值处理。本发明方法结合了U‑net网络与多尺度的匹配滤波算法，与单纯的U‑net网络相比，能够提取到更多的血管特征，具有更高的分割准确性和灵敏度，缓解视网膜血管图像细小血管分割不足以及错误分割的问题。

Description

基于U-net网络的多通道视网膜血管图像分割方法

技术领域

本发明涉及一种图像分割方法，具体是将深度学习中的U-net网络与图像匹配滤波方法相结合进行改进的一种图像分割方法。该方法实际应用在了视网膜血管图像的分割中。

背景技术

目前图像分割的算法主要可划分为两类：监督学习方法和非监督学习方法。非监督学习方法主要是利用图像本身的某些特性，人工设置特征提取的过程与方法，对某些图像的特有属性具有不错的分割结果，但该方法的使用很大程度上依赖于前人经验，在面对不同的图像进行分割处理时，并不具有很好的泛化性。监督学习方法的最主要特征便是需要图像数据集的人工分割结果作为训练的标签，通常该种方法的分割结果具有较好的准确率。

在监督学习算法中，应用于图像分割较为广泛的方法是全卷积神经网络(FCN)。FCN是一个端到端的神经网络，即输出图像尺寸与输入一致，实现了像素级别的分类。其多层卷积层能够在模型训练中提取到不同的特征，而且不同层次的卷积层还可以学习到更为抽象的特征。因此相比于非监督学习方法，该方法能提取到更多的图像特征，能够得到更高的分割准确率。但是传统的FCN存在的问题也是很明显：其网络结构针对像素进行分类，因此不可避免地忽略了像素之间的联系，而且卷积操作的感受野过小，不能很好地利用上下文信息，造成特征提取不完全。

为了克服两种不同分割方法的各自缺点，可以尝试对两种方法进行取长补短。在使用神经网络模型进行训练之前，先对输入图像进行非监督学习方法进行一个预处理，将图像中的全局特征有针对性地提取出来，同时消除掉一些背景噪声，便于后续神经网络模型的特征学习。通过多个通道的不同特征提取，在不改变模型自身结构，不增加训练参数的基础上，可以提取到更多的全局特征，弥补了该特征提取难的缺陷。

发明内容

本发明针对现在图像分割方法的不足之处，将非监督学习方法中的匹配滤波算法用于图像处理中，利用多尺度的特性，在不同通道中提取血管图像中不同尺寸的血管特征；提出一种新的形态学算法用于对血管图像中的伪影等背景噪声进行去除处理；利用基于FCN改进的U-net模型进行多通道特征学习，以获取额外的全局特征，尽可能提高分割的准确率。

按照本发明提供的技术方案，提出了一种基于U-net网络的多通道视网膜血管图像分割方法，包括以下步骤：

步骤1、对训练集图像进行水平翻转、垂直翻转和180度旋转，使数据量扩增4倍；

步骤2、对彩色图像进行图像的预处理操作；

步骤3、构建通道1，对图像进行形态学闭操作用以消除背景噪声；

步骤4、构建通道2，使用大尺度匹配滤波方法和改进的黑帽变换算法处理图像；

步骤5、构建通道3，使用小尺度匹配滤波方法和改进的黑帽变换算法处理图像；

步骤6、将三个通道的图像输入改进的U-net网络进行训练，得到所需的分割网络；

步骤7、融合三个通道的输出分割图，经过后处理后得到最终分割结果。

所述的步骤3中，用于消除背景噪声的一系列操作算法表达式为：

其中，g(x,y)为中间变量，f(x,y)为处理后的图像；I_close(x,y)为进行形态学闭操作后的图像；I(x，y)为步骤2处理后的图；m，n为图像像素尺寸。

所述步骤4、5中的匹配滤波算法表达式如下：

定义二维高斯核函数为：

其中s是高斯核的宽度，l为高斯核的长度，旋转高斯核用以对多方向血管进行滤波。p(x，y)为核函数中一个像素点，旋转矩阵为：

θ_i(0≤θ_i≤p)为第i个核函数的角度，p(x，y)旋转后的坐标值为则第i个模板核函数为：

其中N为模板邻域，取值范围为：N＝{(u，v)，|u|≤3s，|v|≤l/2}。

在实际的算法运用中，常常需要考虑模板滤波器的相关系数均值，如式所示：

其中，A表示模板领域内点的数量，m_i表示第i个滤波器模板中系数的平均值。因此最终模板核函数如式：

所述步骤4、5中的改进黑帽算法表达式如下：

B_hat(I)＝(I(x，y)·b(u，v))-I(x，y)

g(x，y)＝255-I(x，y)-2*B_hat(I)

式中·为形态学闭操作；B_hat(f)为黑帽变换；I(x，y)为原始图像；g(x，y)为最终处理图像。

所述步骤6中，改进的U-net网络结构如下：

步骤6-1：构建编码块。在网络的编码结构，使用了四层编码块，其中每层编码块包括两个卷积层，一个批归一化层(BN)、LReLU激活层以及最大池化层(MaxPooling)。

步骤6-2：构建解码块。在网络的解码结构，与编码结构相对应，使用了四层解码块，其中每层解码块包括两个卷积层，一个上采样层、BN层以及一个跳跃连接层(skipconnection)。

步骤6-3：设置损失函数。为了解决数据集像素分布不平衡的问题，提出一种新的加权损失函数，即将Dice系数与二分类交叉熵函数结合，表达式如下：

Loss＝L_dice+λL_r

L_r为交叉熵函数，表达式如下：

其中TP和TN分别是真阳性和真阴性像素的个数；N_p和N_n分别是分割像素和非分割像素的个数；y是标签值(y＝1，分割目标，y＝0，背景)；p是像素的预测概率值。

L_dice为Dice系数表达式，如下：

其中N是像素的数量值；p(k，i)∈[0，1]，q(k，i)∈[0，1]分别是像素点k类的预测概率和真实标签。

步骤6-4：输出层自适应阈值处理。为了尽可能避免把噪声错分为目标像素的情况出现，采用区域自适应阈值分割算法进行最终处理，表达式如下：

式中b为参数值；(2ω+1)×(2δ+1)为窗口区域面积；T为该区域的阈值；G(x，y)为区域像素点。

所述步骤7的后处理，为了解决存在的孤立噪声点问题。该处理对图像全局进行连通区域面积计算，对连通区域小于25像素面积的区域进行消除处理，即将该处像素重新分类为背景像素。

本发明与现有技术相比，有益效果如下：

利用形态学算法，尽可能消除了视网膜血管图像的背景噪声，有利于神经网络的特征学习；通过多尺度的匹配滤波处理，增强了不同尺寸血管的特征，便于特征提取；利用多通道的特征提取结构，能够提取到更多血管特征，提高网络分割的准确性。与传统U-net网络相比，本发明方法能够分割出更多细小血管，缓解图像分割不足和误分割问题，具有更高的分割灵敏性和准确率。

附图说明

图1为本发明的流程图；

图2为本发明的网络模型结构图；

图3为通道2中的处理结果图；

图4为通道3中的处理结果图；

图5为后处理前后对比图；

图6为各个通道与通道融合后的结果图。

具体实施方式

下面结合具体实施例对本发明做进一步说明。以下描述仅作为示范和接收，并不对本发明作任何形式上的限制。

如图1所示，本发明的实施例实现步骤如下：

步骤1、对现有的公开数据集DRIVE、STARE和CHASE_DB1的训练集进行数据扩增操作，具体为对图像进行水平翻转、垂直翻转和180度旋转，将数据量扩增到原来的4倍。其中，STARE数据集随机选取15张图像为训练集，CHASE_DB1数据集选取前20张图像为训练集。

步骤2、对图像的预处理流程包括：

步骤2-1：对彩色图像进行通道分离，选取血管清晰度较好的绿色通道作为输入图像进行处理，图像尺寸为576×576；

步骤2-2：对步骤2-1获取的灰度图进行限制性直方图均衡化(CLAHE)处理，增强图像对比度，再接着进行Gamma变换(参数值设置为0.6)，调整图像亮度。

步骤3、图像背景去除的一系列操作如下式：

步骤4、对步骤2中获得的图像进行处理，大尺度匹配滤波算法的参数设置：s为1.9，l为10.8,方向数为12。对大尺度血管进行增强后再使用改进的黑帽算法进行背景消除。图3(e)展示了从三个数据集中随机挑选的一张图片的处理结果，可以发现灰度图的伪影噪声得到去除，血管特征得到增强。

步骤5、处理流程同步骤4，小尺度匹配滤波算法参数设置：s为0.1，l为5,方向数为18。处理结果图如图4(e)所示。

步骤6、如图2所示，改进的U-net网络模型具体结构参数如表1所示：

表1：U-net网络结构参数

网络训练采用随机梯度下降(SGD)算法对损失函数进行优化，初始学习率为0.01，epoch设置为40。学习率采用三段式更新的方式，如下式：

步骤7、获得三个通道的分割图后，进行融合处理。为了解决存在的孤立噪声点问题。对图像全局进行连通区域面积计算，随后对连通区域小于25像素面积的区域进行消除处理，即将该处像素重新分类为背景像素。图5给出了后处理的对比图，可以发现红圈中的孤立噪声点得到了有效去除。

在图6中展示了三个数据集图片在三条特征提取通道的最终分割图以及各自的融合图。可以看出各个通道的部分断续血管在融合后重新连接了起来，而且融合图的细小血管数量明显要高于各个单独的通道分割图，血管分割不足的问题得到了一定程度的缓解。为了更直观了解本发明的可行性和优越性，在表2中给出了DRIVE数据集在三个通道的各自性能指标，在表3中对比了专业医师的手工分割指标。

表1:DRIVE数据集的分割结果指标

表3：DRIVE、STARE和CHASE_DB1数据集的分割结果指标

从表2可以看出：通过将三个通道结果进行融合，分割效果明显有了很大的提升，尤其是在灵敏度(Se)指标上，这也从侧面反映了得益于多通道结构，本发明能够分割出更多的细小血管。表3列出了三个数据集第二位专业医师手工分割标签的相应指标，通过对比可以发现，本发明的分割结果在各个指标上基本上要优于手工分割的结果，这也证明了本项发明的可靠性与优越性。

Claims (4)

1.基于U-net网络的多通道视网膜血管图像分割方法，其特征在于包括如下步骤：

步骤1、对训练集图像进行数据量扩增；

步骤2、对彩色图像进行图像的预处理操作；

步骤3、构建通道1，对图像进行形态学闭操作用以消除背景噪声；

步骤4、构建通道2，使用大尺度匹配滤波方法和改进的黑帽变换算法处理图像；

步骤5、构建通道3，使用小尺度匹配滤波方法和改进的黑帽变换算法处理图像；

所述的改进的黑帽变换算法具体为：

B_hat(I)＝(I(x,y)·b(u,v))-I(x,y)

F(x,y)＝255-I(x,y)-2*B_hat(I)

式中·为形态学闭操作；B_hat(I)为黑帽变换；I(x,y)为步骤2处理后的图；F(x,y)为最终处理图像，b(u,v)表示滤波器模板；

步骤6、将三个通道的图像输入改进的U-net网络进行训练，得到所需的分割网络；

所述的步骤6中，构建的U-net网络结构包括依次级联的：

DoubleConv：Conv(ksize＝3,pad＝1)、BN+LreLU、Conv(ksize＝3,pad＝1)、BN+LreLU；

Input；

Encoder block_1：DoubleConv_1、Maxpooling(ksize＝2)；

Encoder block_2：DoubleConv_2、Maxpooling(ksize＝2)；

Encoder block_3：DoubleConv_3、Maxpooling(ksize＝2)；

Encoder block_4：DoubleConv_4、Maxpooling(ksize＝2)；

Encoder block_5：DoubleConv_5、Maxpooling(ksize＝2)；

Decoder block_1：Upsampling(bilinear)ConcatDoubleConv_6；

Decoder block_2：Upsampling(bilinear)ConcatDoubleConv_7；

Decoder block_3：Upsampling(bilinear)ConcatDoubleConv_8；

Decoder block_4：Upsampling(bilinear)ConcatDoubleConv_9；

Output：Conv(ksize＝1)；

在U-net神经网络的训练中，采用新的加权损失函数，表达式如下：

Loss＝L_dice+λL_r

L_r为交叉熵函数，表达式如下：

其中TP和TN分别是真阳性和真阴性像素的个数；N_p和N_n分别是分割像素和非分割像素的个数；y是标签值，其中y＝1表示分割目标，y＝0表示背景；p是像素的预测概率值；

L_dice为Dice系数表达式，如下：

其中N是像素的数量值；p(k,i)∈[0,1]，q(k,i)∈[0,1]分别是像素点k类的预测概率和真实标签；

U-net神经网络在输出层上采用区域自适应阈值分割算法进行最终二值化处理，表达式如下：

式中b为参数值；(2ω+1)×(2δ+1)为窗口区域面积；T为该区域的阈值；G(x,y)为区域像素点；

步骤7、融合三个通道的输出分割图，经过后处理后得到最终分割结果。

2.根据权利要求1所述的基于U-net网络的多通道视网膜血管图像分割方法，其特征在于：

所述的步骤3中，用于消除背景噪声的一系列操作算法表达式为：

其中，g(x,y)为中间变量，f(x,y)为处理后的图像；I_close(x,y)为进行形态学闭操作后的图像；I(x,y)为步骤2处理后的图；m,n为图像像素尺寸。

3.根据权利要求1所述的基于U-net网络的多通道视网膜血管图像分割方法，其特征在于：步骤1对训练集图像进行数据量扩增；扩增方式为对图像进行水平翻转、垂直翻转和多角度度旋转。

4.根据权利要求1所述的基于U-net网络的多通道视网膜血管图像分割方法，其特征在于：步骤4、5中的匹配滤波算法表达式如下：

定义二维高斯核函数为：

其中s是高斯核的宽度，l为高斯核的长度，旋转高斯核用以对多方向血管进行滤波；p(x,y)为核函数中的一个像素点，旋转矩阵为：

θ_i为第i个核函数的角度，0≤θ_i≤p，p为设定的角度值，p(x,y)旋转后的坐标值为 则第i个滤波器核函数为：

其中N为模板邻域，取值范围为：N＝{(u,v),|u|≤3s,|v|≤l/2}；

在实际的算法运用中，常常需要考虑模板滤波器的相关系数均值，如式所示：

其中，A表示模板领域内像素点的数量，m_i表示第i个滤波器模板中系数的平均值；因此最终模板核函数如式：