论文复现之医学图像应用：视网膜血管分割

2019 年 10 月 6 日
笔记

论文复现之医学图像应用：视网膜血管分割

0.导语

今日研究为继续上次论文中的一个内容：U-Net网络，于是找了一篇经典论文，并学习论文及代码解读。在学习U-Net网络后，使用U-Net神经网络提取视网膜纹理血管。

1.论文阅读

论文题目：U-Net Convolutional Networks for Biomedical Image Segmentation。

中文翻译为：用于生物医学图像分割的U-Net卷积网络。

1.1 摘要

之前，在训练一个深度网络需要大量已标注的训练样本。在这篇论文中，提出了一种网络和训练策略。为了更有效的使用标注数据，使用了数据扩张的方法。本片论文的网络分为两部分：收缩路径(contracting path)和扩张路径(expanding path)。

收缩路径主要用于捕捉上下文特征
扩展路径用于定位

该网络获得了赢得了ISBI cell tracking challenge 2015。不仅如此，这个网络非常的快，对一个512×512的图像，使用一块GPU只需要不到一秒的时间。

1.2 介绍

卷积神经网络的典型用途是分类任务，其中输出到图像是单个类别标签。然而，在许多视觉任务中，尤其是在生物医学图像处理中，期望的输出应该包括定位，即：应该将类别标签分配给每个像素。此外，在生物医学任务中，千量级的训练图像通常难以训练。因此，Ciresan等人在滑动窗口设置中训练网络，通过提供围绕该像素的局部区域（补丁）作为输入来预测每个像素的类别标签。

首先，这个网络可以局部化。
其次，补丁方面的训练数据远大于训练图像的数量。

由此产生的网络大幅度赢得了ISBI 2012的EM segmentation challenge。

显然，Ciresan等人的策略有两个缺点。

首先，它非常慢。

因为网络必须分别为每个补丁运行，并且由于补丁重叠而导致大量冗余。

其次，需要在局部准确性和获取整体上下文信息之间取舍。

较大的补丁需要更多的最大池化层来降低局部准确性，而较小的补丁则使网络只能看到很少的上下文。最近的方法提出了一个分类器输出，它考虑了来自多个层的特征。 良好的本地化和上下文的使用是可能的同时。

在这篇文章中，就建立了一个更加优雅的框架，通常被称为“全卷积网络”（FCN)。随后修改并拓展了这个框架，使其可以仅使用少量训练图片就可以工作，获得更高的分割准确率。网络架构如下图所示：

对FCN中的重要修改是：

上采样部分有大量的特征通道，它们允许网络将上下文信息传播到更高分辨率的层。

使得网络架构中扩张路径与收缩路径基本对称，并产生u形结构。

没有任何完全连接的层，分割图仅包含像素，对于该像素，输入图像中的完整上下文是可用的。

该策略允许通过重叠拼贴策略对任意大的图像进行无缝分割。 为了预测图像的边界区域中的像素，通过镜像输入图像来推断丢失的上下文。这种平铺策略对于将网络应用于大图像很重要，否则分辨率将受到GPU内存的限制。

在上述谈了，使用少量数据做训练，在这篇论文中也使用了数据增强方法，对于这篇论文中的任务，论文中通过对可用的训练图像应用弹性变形来使用过多的数据增强。

这使得网络能够学习这种变形的不变性，而不需要在注释图像语料库中看到这些变换。这在生物医学分割中特别重要，因为变形曾经是组织中最常见的变异，并且可以有效地模拟真实变形。 Dosovitskiy等人证明了学习不变性的数据增强在无监督特征学习的范围内的价值。

许多细胞分割任务中的另一个挑战是分离相同类别的触摸物体。

为此，论文中提出使用加权损失：触摸单元之间的分离背景标签在损失函数中获得大的权重。

1.3 网络架构

网络架构如上图所示：它由一条收缩路径（左侧）和一条扩展路径（右侧）组成。

【收缩路径】

收缩路径是典型的卷积网络架构。它的架构是一种重复结构，每次重复中都有2个卷积层和一个pooling层，卷积层中卷积核大小均为3×3，激活函数使用ReLU，两个卷积层之后是一个2×2的步长为2的最大池化层。每一次下采样后我们都把特征通道的数量加倍。

【扩张路径】

扩展路径中的每一步包括对特征映射先进行上采样，然后进行2×2卷积（“上卷积”）。该特征映射将特征通道的数量减半，与收缩路径中相应裁剪的特征拼接起来，以及两个3×3卷积，每个卷积都有一个ReLU。由于每个卷积中边界像素的丢失，裁剪是必要的。在最后一层，使用1×1卷积来将每个64通道特征图映射到特定深度的结果。

【卷积层总数】

网络总共有23个卷积层。

1.4 训练

论文中采用随机梯度下降法(SGD)训练。由于没有使用0填补的卷积，输出图像比输入小一个恒定的边界宽度。为了最大限度地降低开销并最大限度地利用GPU内存，我们倾向于在较大批量的情况下使用较大的输入切片，从而将批量减少为单个图像。因此，我们使用高动量（0.99），以便大量先前看到的训练样本确定当前优化步骤中的更新。

损失函数使用逐像素的 softmax 函数和交叉熵损失函数的结合
Softmax函数：

a_k(x)表示在feature maps中的的channel=k的feature map像素位置为x的激活值。K表示总类别数，pk(x)表示似然函数。

交叉熵函数：

其中l是每个像素的真实标签，w是权重地图，表示训练中某些像素更加重要。

为了使某些像素点更加重要，我们在公式中引入了w(x)。我们对每一张标注图像预计算了一个权重图，来补偿训练集中每类像素的不同频率，使网络更注重学习相互接触的细胞之间的小的分割边界。我们使用形态学操作计算分割边界。权重图计算公式如下：

其中w_c是权重地图来平衡类像素的频率。d1表示最近单元边界的距离。d2表示到第二最近单元的边界的距离。文中设置w_0=10,σ≈5pixels。

训练需要标注对应的mask,就是类别的区域标记。

网络中权重的初始化：在具有许多卷积层和通过网络的不同路径的深层网络中，权重的良好初始化非常重要。否则，网络的某些部分可能会过度激活，而其他部分则无法提供。理想情况下，应调整初始权值，使网络中的每个特征映射具有近似单位差异。对于具有我们的体系结构（交替卷积和ReLU层）的网络，这可以通过从具有标准偏差(2/N)^0.5的高斯分布中绘制初始权重来实现，其中N为每个神经元的输入节点数量。例如3×3的64通道的卷积层N=3x3x64=576。

1.5 数据增加

在只有少量样本的情况况下，要想尽可能的让网络获得不变性和鲁棒性，数据增强是必不可少的。在显微图像的情况下，我们主要需要移位和旋转不变性以及对变形和灰度值变化的鲁棒性。论文中使用随机位移矢量在粗糙的3×3网格上产生平滑形变。位移是从10像素标准偏差的高斯分布中采样的。然后使用双三次插值计算每个像素的位移。在压缩路径末端的退出层执行进一步的隐式数据增强。

最后论文的实验部分，这里直接在DRIVE数据库上做实验！

2.视网膜血管分割实验

2.1 实验任务

实验任务：使用U-Net神经网络提取纹理血管。

为什么要做这个，有什么实际意义？

临床实验中我们要能够更好的对眼部血管等进行检测、分类等操作，我们首先要做的就是对眼底图像中的血管进行分割，保证最大限度的分割出眼部的血管。从而方便后续对血管部分的操作。

2.2 数据集

实验数据集为：DRIVE数据集。

数据集下载地址： http://www.isi.uu.nl/Research/Databases/DRIVE/

【数据集识别】

数据集分为训练集与测试集，而每个都分为如下目录：

1st_manual目录下的数据集为：手工分好的血管图像，如下图所示：

images目录下的数据集为眼底图像，如下图所示：

mask目录下的数据集为眼部轮廓图像，如下图所示：

从上面可以看出数据集非常少，总共只有20幅图片，为此我们需要做预处理增大数据集。

【U-Net网络优势】

在上述U-net论文中提到U-Net网络可以针对很少的数据集来进行语义分割，比如我们这个眼球血管分割就是用了20张图片来训练就可以达到很好的效果。而且我们这种眼球血管，或者指静脉，指纹之类的提取特征或者血管静脉在U-net网络里就是一个二分类问题。而本文用的U-net网络来实现这个二分类就只需要二十张图片来作为数据集。

2.2 运行环境及库

复现论文代码：

官方：https://github.com/orobix/retina-unet 非官方：https://github.com/CNU105/Retina-Unet

numpy >= 1.11.1
keras >= 2.1.0
PIL >=1.1.7
opencv >=2.4.10
h5py >=2.6.0
configparser >=3.5.0b2
scikit-learn >= 0.17.1

上面是官方给出的运行库环境。实际上当你安装上后还是有问题的！

还需要安装如下内容：

graphviz
pydot
pydot-ng

如果graphviz安装报错，则需要：

sudo apt-get install graphviz或者yum安装

但是我的是服务器上，没有root权限，怎么安装呢？

这就需要源码编译，在linux系统中，源码编译基本可以解决一切安装问题！

由于没有安装graphviz，同时需要root权限，所以安装就出现了问题，然后只能通过源码编译，结果源码编译最新版也出现了问题，说是最新版的bug，我就转向系统的yum本地用户安装：

# 查看yum中是否有你需要的包  yum list 'graphviz*'  # 下载rpm包(上述输出：graphviz.x86_64)  yumdownloader graphviz.x86_64  # 解压RPM包  rpm2cpio graphviz-2.30.1-19.el7.x86_64.rpm |cpio -idvm  # 添加环境变量  vim  ~/.bashrc  export PATH=$PATH:$HOME/usr/bin/  source ~/.bashrc  # 备注：这种安装会在当前目录下新建一个usr目录，安装的软件全在usr下的bin目录下！

安装好上述就可以了。

但是我的服务器上没python，我也没root权限，如何安装python呢?

那就是上述说的源码编译，指定源码编译路径，然后添加环境变量即可！

过程中重要两点：

（1）编译路径指定

./configure --prefix=/home/city/python3/

编译后，直接make，然后make install即可。

（2）配置环境变量

vim  ~/.bashrc  export PATH=$PATH:$HOME/home/city/python3/  source ~/.bashrc

至此，通过源码编译与本地用户安装解决了环境问题！

2.3 开始实验

【数据预处理】

上面说数据集的时候，提到数据集非常少，那么在实验前就得做数据预处理，而该数据预处理做的工作就是：读取原来的训练集与测试集，重新生成数据文件，并保存为hdf5文件。

对应到代码中就是运行下面代码：

python3 prepare_datasets_DRIVE.py

运行该代码后，会生成如下目录及文件：

【配置文件】

在代码中有个configuration.txt配置文件，该配置文件主要用来配置数据集目录，训练及测试设置！

【开始训练】

运行代码

在官网中是让我们运行下面代码：

python3 run_training.py

而实际上直接可以运行：

python3 ./src/retinaNN_training.py

因为，上述代码实际上就是获取我们的实验设置参数，判断是否是GPU训练而已。如果自己会后台运行，直接敲命令即可！

通过运行retinaNN_training.py我们的实验便开始进入训练过程了，训练过程非常漫长，本次运行GPU共8块，运行时长达到6个h。

模型查看

网络定义每行这个小括号填inputs是代表这层模型连接在inputs之后

inputs = Input(shape=(n_ch,patch_height,patch_width))  # 2D 卷积层  conv1 = Conv2D(32, (3, 3), activation='relu', padding='same',data_format='channels_first')(inputs)

完整模型注释及代码：

def get_unet(n_ch,patch_height,patch_width):      # 第一层      inputs = Input(shape=(n_ch,patch_height,patch_width))      # 2D 卷积层      conv1 = Conv2D(32, (3, 3), activation='relu', padding='same',data_format='channels_first')(inputs)      # Dropout抛弃一些参数防止过拟合      conv1 = Dropout(0.2)(conv1)      conv1 = Conv2D(32, (3, 3), activation='relu', padding='same',data_format='channels_first')(conv1)      pool1 = MaxPooling2D((2, 2))(conv1)        # 第二层      conv2 = Conv2D(64, (3, 3), activation='relu', padding='same',data_format='channels_first')(pool1)      conv2 = Dropout(0.2)(conv2)      conv2 = Conv2D(64, (3, 3), activation='relu', padding='same',data_format='channels_first')(conv2)      pool2 = MaxPooling2D((2, 2))(conv2)        # 第三层      conv3 = Conv2D(128, (3, 3), activation='relu', padding='same',data_format='channels_first')(pool2)      conv3 = Dropout(0.2)(conv3)      conv3 = Conv2D(128, (3, 3), activation='relu', padding='same',data_format='channels_first')(conv3)      # 上采样函数，size(2,2)其实就等于将原图放大四倍（水平两倍，垂直两倍）例如32*32 变成 64*64的图像      up1 = UpSampling2D(size=(2, 2))(conv3)        # 两层拼接      up1 = concatenate([conv2,up1],axis=1)        conv4 = Conv2D(64, (3, 3), activation='relu', padding='same',data_format='channels_first')(up1)      conv4 = Dropout(0.2)(conv4)      conv4 = Conv2D(64, (3, 3), activation='relu', padding='same',data_format='channels_first')(conv4)      up2 = UpSampling2D(size=(2, 2))(conv4)        # # 两层拼接      up2 = concatenate([conv1,up2], axis=1)        conv5 = Conv2D(32, (3, 3), activation='relu', padding='same',data_format='channels_first')(up2)      conv5 = Dropout(0.2)(conv5)      conv5 = Conv2D(32, (3, 3), activation='relu', padding='same',data_format='channels_first')(conv5)      conv6 = Conv2D(2, (1, 1), activation='relu',padding='same',data_format='channels_first')(conv5)      conv6 = core.Reshape((2,patch_height*patch_width))(conv6)      # Permute层将输入的维度按照给定模式进行重排      conv6 = core.Permute((2,1))(conv6)      conv7 = core.Activation('softmax')(conv6)        # 将模型的输入和输出给Model函数就会自己组建模型运行图结构      model = Model(inputs=inputs, outputs=conv7)      # sgd = SGD(lr=0.01, decay=1e-6, momentum=0.3, nesterov=False)      # 配置模型      model.compile(optimizer='sgd', loss='categorical_crossentropy',metrics=['accuracy'])      return model

模型训练代码流程：

首先定义模型
读取配置文件
加载数据

这一步非常重要，对读入内存准备开始训练的图像数据进行一些增强之类的处理。这里对其进行了，直方图均衡化，数据标准化，并且压缩像素值到0-1，将其的一个数据符合标准正态分布。(这个处理代码在lib目录下的extract_patches.py中)。

前面我们知道数据集非常少，那怎么让训练的数据集增大呢？

这里就是做这个处理，即对每张图片取patch时，除了正常的每个patch每个patch移动的取之外，我们还在数据范围内进行随机取patch，这样虽然各个patch之间会有一部分数据是相同的，但是这对于网络而言，你传入的也是一个新的东西，网络能从中提取到的特征也更多了。这一步的目的其实就是在有限的数据集中进行一些数据扩充，这也是在神经网络训练中常用的手段了。

存储随机组合的patch

左图为随机原图，右图为mask图

存储模型

【测试及评估】

python3 src/retinaNN_predict.py

由于服务器上没有可视化页面，所以当调用matlibplot时候，会报如下错误，那么就需要后台运行，不显示！

no display name and no $DISPLAY environment variable

解决办法是：

加入下面两行代码即可！

import matplotlib  matplotlib.use('Agg')

在本代码中，做了如下工作：

同训练集一样取patch
利用训练好的模型加载测试数据
存储并可视化结果
评估结果

从左到右，依次为真实，原图，预测。

从上到下，依次为真实，原图，预测。