CN109272507A - 基于结构随机森林模型的相干光断层图像的层分割方法 - Google Patents

Abstract

本发明公开了一种基于结构随机森林模型的相干光断层图像的层分割方法，包括步骤：S1、图像预处理：对输入图像进行像素值转化和去噪处理；S2、特征提取：先提取手工特征，再训练一个残差网络提取深度特征，然后整合手工特征和深度特征；S3、通过结构随机森林获取层轮廓的概率图：使用训练图像集中数据训练能预测层轮廓的结构随机森林，然后预测测试图像中的层轮廓的概率图；S4、利用最短路径进行层分割：利用最短路径获得具体的层边界位置。该方法结合深度特征和手工特征的优点，手工特征和深度特征的组合可以训练一个更高效的随机结构森林，使得分割性能得到加强；该方法可以分割正常的视网膜图像，也可以分割病变的视网膜图像，具有较高的灵活性和适应性。

Description

基于结构随机森林模型的相干光断层图像的层分割方法

技术领域

本发明涉及相干光断层图像处理技术领域，特别是涉及一种基于结构随机森林模型的相干光断层图像的层分割方法。

背景技术

光学相干断层扫描成像技术(Optical Coherence Tomography,OCT)在眼科中是一种重要的非侵入性的成像方式。时域OCT(TD-OCT)是第一代OCT技术，每秒可获得400张A扫描图，轴向分辨率范围为8到10微米。当前商用的OCT 技术是谱域OCT(SD-OCT)，成像速度比时域OCT快10到100倍，轴向分辨率范围为3到7微米，可以获取视网膜和角膜的切面信息，由于成像速度快，分辨率高，谱域OCT早期用于疾病的诊断和预测。

如图6所示，视网膜的层结构依次为：Vitreous(玻璃体)、NFL(神经纤维层)、GCL+IPL(神经节细胞层和内丛状层)、INL(内核层)、OPL(外丛状层)、 ONL-ISM(外核层和内段骨髓)、ISE(内段椭球)、OS+RPE(光感受器外节和视网膜色素上皮)和Choroid(脉络膜)。层分割是眼科疾病诊断的关键步骤，许多眼科疾病会导致视网膜层在解剖学上的变化，如年龄相关性黄斑病变(AMD)，斯特格病变(SD)，视网膜色素病变(RP)和糖尿病黄斑水肿(DME)，这也使得视网膜层分割非常困难。手工分割已经应用在许多OCT研究中。然而，手工分割是一个耗时的过程，并且具有很强的主观性。为了解决这个问题，相关研究人员已经提出了许多层分割方法，如主动轮廓方法、水平集方法、动态规划方法和图割方法。近来，越来越多的基于机器学习的方法被提出，其中，基于动态规划和基于机器学习的方法已经得到广泛应用。Chiu提出了一种基于最短路径的自动层分割方法，Karri针对视网膜图像中特定层边缘检测提出一种结构化学习算法，增强了传统的基于图的层分割效果。然而，所有这些方法仅利用手工提取的特征，对于许多情况不适合，例如形变与噪声。

近年来，深度学习(DL)在许多低层次和高层次的计算机视觉(CV)任务中都取得了巨大的成功，如图像去噪、图像分类、图像分割、目标检测与识别，也被用于医学图像处理。深度特征取代手工特征(像HOG与SIFT)的趋势已经出现，于是将深度学习技术引入到OCT图像的层分割是合理且可行的。Fang通过卷积神经网络(CNN)将图搜索算法应用到概率输出中，从而提出一种视网膜的九层分割方法。ReLayNet是一种类似完全卷积网络的U-net，可以用于点对点训练中的层分割。研究表明手工提取特征和深度特征是互补的，两者结合在一起通常得到更好的结果，其中网络的深度是一个关键的参数，随着梯度的增加，会出现梯度消失和爆炸的问题，妨碍收敛。He提出一种残差学习框架加快网络的训练，实质上比其他网络更深。

通常，OCT B扫描图包含噪声，不同层之间的区分度不高；并且，形变区域会使得手工特征提取困难。

发明内容

本发明为了克服以上技术的不足，提供了一种基于结构随机森林模型的相干光断层图像的层分割方法。

术语解释：

1、SGD：Stochastic Gradient Descent，指随机梯度下降算法。

2、BM3D算法：三维块匹配算法。

本发明克服其技术问题所采用的技术方案是：

一种基于结构随机森林模型的相干光断层图像的层分割方法，包括如下步骤：

S1、图像预处理：

对输入图像进行像素值转化和去噪处理；

S2、特征提取：

先提取手工特征，再训练一个残差网络提取深度特征，然后整合手工特征和深度特征；

S3、通过结构随机森林获取层轮廓的概率图：

使用训练图像集中数据训练能预测层轮廓的结构随机森林，然后预测测试图像中的层轮廓的概率图；

S4、利用最短路径进行层分割：

在步骤S3获得的概率图上，利用最短路径获得具体的层边界位置。

本发明优选的，所述步骤S1中，将输入图像的像素值范围从[0，255]转化到[0，1]，采用BM3D算法对所有图像进行去噪处理。

本发明优选的，所述手工特征包括四个方向的梯度特征、方向梯度直方图特征和Gabor特征；提取手工特征的具体步骤如下：

采用高斯滤波器的一阶和二阶导数生成高斯特征，得到五个高斯特征图，包括两个一阶高斯特征图和三个二阶高斯特征图；

采用Gabor滤波器提取了在180°范围内的最大响应值作为Gabor特征，得到一个Gabor特征图；

输入图像采用双线性向下采样插值来模拟多尺度，然后在降采样后的图像上计算梯度幅值图、Gabor特征图和五个高斯特征图；

所有手工特征图总共是15维特征图，包含一个灰度图和在两个尺度上提取的特征图，所述每个尺度特征包括一个梯度幅度特征图、一个Gabor特征图、两个一阶高斯特征图和三个二阶高斯特征图，将所有特征图沿着第三个维度连接起来，形成最终的手工特征。

本发明优选的，步骤S2中，训练残差网络提取深度特征时，从数据集的每个图像块中采样得到训练数据和测试数据，用于训练卷积神经网络CNN模型；在训练残差网络模型CNN时，采用SGD算法迭代最小化损失，用于优化的交叉熵损失函数被定义如下：

其中Indicator(·)是指示函数，N是样本的数量，K＝8是样本类数，lⁱ是第i个样本的真实标签，是CNN模型的输出概率，第i个样本的标签是j；

对于一个m×n大小的OCT图像，提取m×n个中心在每个像素的块，对于每个块，将最后一层softmax层作为深度特征。

本发明优选的，步骤S3中，训练结构随机森林的层轮廓时，使用一个输入图像、与所述输入图像对应的标签图像以及与所述输入图像对应的轮廓图像，所述标签图像是通过将输入图像中的每个像素替换为对应的层标签构建而成，所述轮廓图像是通过替换输入图像中的上边界像素为对应的层标签，同时抑制上边界像素以外的像素而成；对于输入图像和标签图像的每一层，执行指导抽样，针对轮廓训练结构随机森林，依此步骤重复每一层的训练，层数用c表示，其中c＝{1,2，……8}，保存每个层轮廓的结构随机森林。

本发明优选的，步骤S4中，采用动态规划算法确定从图像左上方到右下方的最短路径，对于分割第c层边界，每一张概率图被视为一张图，图中边的权重是基于相应的概率图计算，公式如下：

其中a是一个像素，b是a的邻接像素，和是像素a和b属于边界l的概率，w_m是最小边缘权重，被用来寻找最短路径。

本发明的有益效果是：

本发明的方法结合深度特征和手工提取的特征来训练结构随机森林。手工特征包括四个方向的梯度特征、方向梯度直方图(HOG)特征和Gabor特征，深度特征来自于深度残差网络，该网络模型利用有标签的数据通过监督学习得到。提取出图像的手工特征和深度特征后，结合这些特征训练一个随机结构森林，利用训练好的随机结构森林，生成每层的轮廓概率图，然后使用最短路径获取最终的分割结果。本发明的方法具有如下优点：

(1)高效性。本方法结合深度特征和手工特征的优点，手工特征和深度特征的组合可以训练一个更高效的随机结构森林，使得分割性能得到加强。

(2)灵活性。本方法可以分割正常的视网膜图像，也可以分割病变的视网膜图像，具有较高的灵活性和适应性。

附图说明

图1是本发明训练过程的流程示意图。

图2为本发明预测过程的流程示意图。

图3为本发明实施例所使用的深度特征提取器模型。

图4为本发明实施例的结构随机森林训练过程。

图5为本发明实施例的结构随机森林预测过程。

图6是视网膜层结构示意图。

图7是本发明实施例选取的一张视网膜图。

图8是本发明实施例专家一标注的结果。

图9是本发明实施例专家二标注的结果。

图10是本发明实施例分割的结果图。

具体实施方式

为了便于本领域人员更好的理解本发明，下面结合实施例对本发明做进一步详细说明，下述仅是示例性的不限定本发明的保护范围。

本发明提供了一种基于结构随机森林模型的相干光断层图像的层分割方法，下面结合附图1-10说明其具体实施过程。

步骤一、数据预处理

本发明使用的图像数据来自Duke DME数据集，包括若干张视网膜OCT B扫描图，每张扫描图由多名专家标注8个层结构的上边界，本实施例中图像是由 10个人提供的110张OCT B扫描图组成，每张扫描图由两名专家标注8个层结构的上边界。

先对输入图像的像素值进行转换：将输入图像(I∈R^m×n)的像素值范围从 [0，255]转化到[0，1]，图像采集协议可能导致图像边界处的值饱和，所有在I 中等于1的值替换为0.01。

然后，对所有图像通过基于块的图像去噪方法处理，本实施例采用BM3D算法进行去噪处理，所述BM3D算法：它首先把图像分成一定大小的块，根据图像块之间的相似性，把具有相似结构的二维图像块组合在一起形成三维数组，然后用联合滤波的方法对这些三维数组进行处理，最后，通过逆变换，把处理后的结果返回到原图像中，从而得到去噪后的图像。

步骤二、特征提取

1、手工特征的提取：

在视网膜OCT B扫描图中，相邻两层之间存在灰度的变化，这是最显著的特征，因此，手工特征包括图像的强度和梯度特征。OCT图像被认为是一幅灰度图。在每个像素上计算的归一化后的梯度幅值可以形成一个幅度图。另外，Gabor 滤波器和高斯滤波器用来生成其他两种特征。

本实施例中，我们使用高斯滤波器的一阶和二阶导数生成高斯特征，这样就会有五个特征图(两个一阶高斯特征图和三个二阶高斯特征图)，sigma参数设置为0.5。

Gabor滤波器是一种属于窗口傅里叶变换的线性滤波，它在提取目标局部空间和频域信息具有良好的性能。本实施例中，我们提取了在180°范围内的最大响应值作为Gabor特征，到一个Gabor特征图。

此外，输入图像采用双线性向下采样插值来模拟多尺度，然后在降采样后的图像上计算梯度幅值图、Gabor特征图和五个高斯特征图。和原始图像灰度值一起总共15维特征图(所有的手工特征图包含一个灰度图和在两个尺度上提取的特征图。每个尺度特征包括一个梯度幅度特征图、一个Gabor特征图、两个一阶高斯特征图和三个二阶高斯特征图)形成。最后，所有特征图沿着第三个维度连接起来，形成最终的手工特征，大小为m×n×15。

2、深度特征的提取：

为了生成深度特征，首先需要训练一个残差网络，残差网络是指包含残差块的卷积神经网络，加入残差块的网络可以有效避免梯度爆炸或消失现象。卷积神经网络CNN模型中用到的训练数据和测试数据是从数据集的每个图像中采样得到。具体来说，从每个训练图像中提取65×65的图像块，每个块的标签取决于该65×65的图像块的中心像素，每个块的标签与该65×65图像块的中心像素的标签相等，边界上图像块中心像素的标签(ILM的上边界到RPE的下边界) 被分配到1-8，边界以外的块标签为0，得到9类样本，所述9类样本中，每一类的样本数量都是相同的。

一个典型的卷积神经网络CNN包括卷积层、池化层、批处理归一化层和ReLU 层，基线的卷积层包括19个卷积层和一个全连接层，第一个卷积层使用16个卷积核大小为3×3，下采样的步幅设置为2，另外18个卷积层中，6个卷积层使用16个卷积核，大小为3×3，6个卷积层使用32个卷积核，大小为3×3，6 个卷积层使用64个卷积核，大小为3×3，所有卷积层都经过批处理归一化层和 ReLU层处理，最后的一个3×3的卷积层经过8×8的池化层处理后，附加9个神经元的全连接层，最后一层是softmax层。

在训练残差网络模型CNN时，采用SGD算法迭代最小化损失，用于优化的交叉熵损失函数被定义如下：

其中Indicator(·)是指示函数，N是样本的数量，K＝8是样本类数，lⁱ是第i个样本的真实标签，是CNN模型的输出概率，第i个样本的标签是j；

特别地，对于一个m×n大小的OCT图像，提取m×n个中心在每个像素的块，块大小为65×65。对于每个块，将最后一层softmax层作为深度特征。具体说来，对于给定的OCT图像中的每个像素，生成一个9维的特征向量，这些特征向量被用来生成一个9维的深度特征，大小为m×n×9。

为了获得更有效的图像特征，手工特征和深度特征需要整合，我们沿着第三维将手工特征和深度特征整合，手工特征大小为m×n×15，深度特征大小为 m×n×9，如此，整合后图像特征大小为m×n×24。

步骤三、通过结构随机森林获取层轮廓的概率图

该部分分为两个阶段，第一阶段是训练结构随机森林的层轮廓，第二阶段是预测测试图像中的每层上边界的层轮廓的概率图。层轮廓预测问题可以被视为多分类问题或二分类问题。本实施例中，我们将它视为二分类问题，所以对于每个层轮廓需要训练结构随机森林。结构随机森林输出八个层轮廓概率图，这将用于基于图的分割方法。

训练结构随机森林的层轮廓时，使用一个输入图像、与所述输入图像对应的标签图像以及与所述输入图像对应的轮廓图像。

首先进行数据预处理，与步骤一的预处理步骤相同；然后将手工特征和深度特征整合成最终的图像特征。所述标签图像是通过将输入图像中的每个像素替换为对应的层标签构建而成，所述轮廓图像是通过替换输入图像中的上边界像素为对应的层标签，同时抑制上边界像素以外的像素而成。对于输入图像和标签图像的每一层，执行指导抽样，针对轮廓训练结构随机森林，依此步骤重复每一层的训练，层数用c表示，其中c＝{1,2，……8}，保存每个层轮廓的结构随机森林。

在层轮廓的概率图预测阶段，也要进行数据预处理和特征提取，数据预处理和特征提取的过程与训练阶段相同。通过训练后的结构随机森林，整合后的特征可以用来预测层轮廓概率图。第n层上边界的预测结果是一张概率图，大小为m×n。

步骤四、利用最短路径进行层分割

给定一个输入图像，大小为m×n，每个训练后的结构随机森林输出相应的概率图，大小为m×n。如果直接确定像素属于某一层边界的概率可能不是最优的，因为它不能保证最终层分割的连通性和光滑性。因此，采用一种改进的基于图的动态规划算法确定从图像左上方到右下方的最短路径，用于精确分割边界。对于分割第c层边界，每一张概率图被视为一张图，图中边的权重是基于相应的概率图计算，公式如下：

其中a是一个像素，b是a的邻接像素，和是像素a和b属于边界l的概率，w_m是最小边缘权重，被用来寻找最短路径。

以上仅描述了本发明的基本原理和优选实施方式，本领域人员可以根据上述描述做出许多变化和改进，这些变化和改进应该属于本发明的保护范围。

Claims (6)

1.一种基于结构随机森林模型的相干光断层图像的层分割方法，其特征在于，包括如下步骤：

S1、图像预处理：

对输入图像进行像素值转化和去噪处理；

S2、特征提取：

先提取手工特征，再训练一个残差网络提取深度特征，然后整合手工特征和深度特征；

S3、通过结构随机森林获取层轮廓的概率图：

使用训练图像集中数据训练能预测层轮廓的结构随机森林，然后预测测试图像中的层轮廓的概率图；

S4、利用最短路径进行层分割：

在步骤S3获得的概率图上，利用最短路径获得具体的层边界位置。

2.根据权利要求1所述的方法，其特征在于，所述步骤S1中，将输入图像的像素值范围从[0，255]转化到[0，1]，采用BM3D算法对所有图像进行去噪处理。

3.根据权利要求1或2所述的方法，其特征在于，步骤S2中，所述手工特征包括四个方向的梯度特征、方向梯度直方图特征和Gabor特征；提取手工特征的具体步骤如下：

采用高斯滤波器的一阶和二阶导数生成高斯特征，得到五个高斯特征图，包括两个一阶高斯特征图和三个二阶高斯特征图；

采用Gabor滤波器提取了在180°范围内的最大响应值作为Gabor特征，得到一个Gabor特征图；

输入图像采用双线性向下采样插值来模拟多尺度，然后在降采样后的图像上计算梯度幅值图、Gabor特征图和五个高斯特征图；

所有手工特征图总共是15维特征图，包含一个灰度图和在两个尺度上提取的特征图，所述每个尺度特征包括一个梯度幅度特征图、一个Gabor特征图、两个一阶高斯特征图和三个二阶高斯特征图，将所有特征图沿着第三个维度连接起来，形成最终的手工特征。

4.根据权利要求3所述的方法，其特征在于，步骤S2中，训练残差网络提取深度特征时，从数据集的每个图像块中采样得到训练数据和测试数据，用于训练卷积神经网络CNN模型；在训练残差网络模型CNN时，采用SGD算法迭代最小化损失，用于优化的交叉熵损失函数被定义如下：

其中Indicator(·)是指示函数，N是样本的数量，K＝8是样本类数，lⁱ是第i个样本的真实标签，是CNN模型的输出概率，第i个样本的标签是j；

对于一个m×n大小的OCT图像，提取m×n个中心在每个像素的块，对于每个块，将最后一层softmax层作为深度特征。

5.根据权利要求1或2或4所述的方法，其特征在于，步骤S3中，训练结构随机森林的层轮廓时，使用一个输入图像、与所述输入图像对应的标签图像以及与所述输入图像对应的轮廓图像，所述标签图像是通过将输入图像中的每个像素替换为对应的层标签构建而成，所述轮廓图像是通过替换输入图像中的上边界像素为对应的层标签，同时抑制上边界像素以外的像素而成；对于输入图像和标签图像的每一层，执行指导抽样，针对轮廓训练结构随机森林，依此步骤重复每一层的训练，层数用c表示，其中c＝{1,2，……8}，保存每个层轮廓的结构随机森林。

6.根据权利要求5所述的方法，其特征在于，步骤S4中，采用动态规划算法确定从图像左上方到右下方的最短路径，对于分割第c层边界，每一张概率图被视为一张图，图中边的权重是基于相应的概率图计算，公式如下：

其中a是一个像素，b是a的邻接像素，和是像素a和b属于边界l的概率，w_m是最小边缘权重，被用来寻找最短路径。