CN107392918B - 基于随机森林与复合活性曲线的oct图像层分割方法 - Google Patents

Abstract

本发明涉及一种基于随机森林与复合活性曲线的OCT图像层分割方法，为了精确的分割视网膜层与积液而设计。本发明基于随机森林与复合活性曲线的OCT图像层分割方法，包括：得到OCT图像特征训练随机森林分类器；复合活性曲线算法获得最终的SF1；提取OCT图像的24个特征，使用随机森林分类器，本发明基于随机森林与复合活性曲线的OCT图像层分割方法，操作简单、检测结果准确。克服现有的对于病变OCT图像分割算法识别率较低、分割效果较差等问题。

Description

基于随机森林与复合活性曲线的OCT图像层分割方法

技术领域

本发明属于医学图像处理算法领域，具体涉及一种基于随机森林与复合活性曲线的OCT图像层分割方法。

背景技术

中心浆液性视网膜病变是严重复杂的视网膜疾病，容易导致失明。中心浆液性视网膜病变发生主要表现为交错区视网膜下浆液的积累，也可能导致视网膜色素上皮脱离。除浆液积聚外，视网膜色素上皮脱离可能出现在浆液下，也可发生在黄斑中心附近。这些液体导致视网膜层肿胀，视网膜层的厚度和光学强度可能会突然改变。因此，中心浆液性视网膜病变是一种常见的黄斑疾病，而黄斑负责中央视觉。因此关于中心浆液性视网膜病变的定量分析研究在视网膜的研究中具有非常重要的意义。

光学相干断层扫描(OCT)是一种无创、非接触的视网膜异常形态成像技术。目前，OCT设备有了很大的进步，使得量化视网膜疾病更准确的量化、异常的视网膜层厚度计算研究更加成为可能。但是，OCT图像中视网膜层与浆液分割技术缺乏、分割精度不高，原因在于：首先，视网膜层的形状复杂的，它的内部边界是非光滑的，有超过十二层结构。其次，有可能包含病变结构，如浆液，视网膜脱离引起的积液。这导致OCT图像视网膜层之间低对比度和模糊边界，也视网膜层结构变化很大。因此，使用传统的表面检测方法如图搜索算法，层分割可能会失败，同时，积液分割使用传统的方法，如区域增长，也可能很容易泄漏到邻域。

鉴于上述，本设计人积极加以研究创新，以期创设一种基于随机森林与复合活性曲线的OCT图像层分割方法，使其更具有产业上的利用价值。

发明内容

为解决上述技术问题，本发明的目的是为提供一种分割准确性高的基于随机森林与复合活性曲线的OCT图像层分割方法

为达到上述发明目的，本发明基于随机森林与复合活性曲线的OCT图像层分割方法，包括：

训练随机森林分类器；在训练随机森林分类器时，将OCT图像的多个层结构分割为了n个分类标记区域，各区域分别具有上表面，各区域的上表面层由上至下依次为SF1、SF2、SF3……SFn；提取了OCT图像中的m个特征；

确定待分割OCT图像上表面SF1的初始形状，使用复合活性曲线算法得到上表面SF1的最终形状；

提取待分割OCT图像与训练随机森林分类器时相同的m个特征，使用随机森林分类器，对待分割OCT图像进行分类标记为n个区域，找到区域2至区域n中N个区域上表面的初始形状，基于各区域上表面的初始形状分别各自使用复合活性曲线算法，得到各区域上表面的最终形状，其中n、N均为正整数，N≤n-1；

当N＜n-1时，还包括：将一区域的上表面的最终形状向上或向下平移预定像素点得到剩余其他一区域的上表面的初始形状，基于该剩余其他一区域的上表面的初始形状，使用复合活性算法得到该剩余其他一区域的上表面的最终形状；

还包括：判断一区域内是否包含OCT图像的至少两层结构，若包含至少两层结构，则首先得到其中一层结构的上表面的初始形状，使用随机森林活性算法得到该层结构的上表面的最终形状，然后该层结构的最终形状向上或向下平移预定像素点得到剩余其他一层结构的上表面的初始形状，基于该剩余其他层结构的上表面的初始形状，使用复合活性算法得到该剩余其他一层结构的上表面的最终形状；

其中，复合活性曲线算法具体包括：

步骤A1：对待分割的OCT图像使用各向异性滤波后的OCT图像

并计算亮层结构的多尺度响应响应

其中，

表示体素的坐标，σ_t表示尺度；

步骤A2：分别从

和

依次提取x-z平面图像，得到各自的二维图像I_2d和L_2d；提取出初始表面SF在当前切面中初始曲线

如果待分割曲面在OCT图像是由暗带到亮带交界面，则上下反转二维图像I_2d和L_2d；

如果待分割曲面在OCT图像是由亮带到暗带交界面，则不反转；

步骤A3：分别对二维图像I_2d和L_2d建立图G^F＝(V^F,E^F)和G^B＝(V^B,E^B)，其中V表示图的结点集合，E表示图的有向边集合；对于二维图像来说，图的结点是四个相邻像素交点，图的边是两个相邻像素p,q相接的边，有向边的右边像素表示内部、左边像素表示外部，按如下公式计算E^F和E^B向边e_ij的权值：

其中，f是特征函数，c_k是传递函数，ω_k是权值系数，n_l是边特征总数；n_k是传递函数总数；复合图G^H中的结点V^H与边E^H和图G^F中的结点V^F与边E^F、G^B中的结点V^B与边E^B一一对应，根据E^F和E^B计算复合图G^H＝(V^H,E^H)边的权值：

其中

表示点积；

将初始曲线

看成更新曲线

能量值E'←∞,E_lw←0，采用如下方式更新：

步骤B1：改变能量值E'←E_lw,E_lw←0，按间距lw_s从曲线

采样得到两个结点a_o和a_o+1，

a_x表示当前二维图像宽度；

步骤B2：在复合图G^H中，使用Dijkstra算法找到a_o到a_o+1的最小代价路径

其中，n₀表示当前路径途径结点总数；计算当前路径的代价值：

并更新代价值E_lw←E_lw+e_lw；代价路径

中等距采样得到lw_s结点更新曲线

步骤B3：更新采样步长lw_s←lw_s/2，如果|E_lw-E'|＞ΔE并且lw_s＞1，则返回步骤1；否则，平滑已更新曲线

如果输入图像上下反转，则把曲线

反转回去；最后，输出已更新曲线

进一步地，训练随机森林分类器的训练包括：

将将中心浆液性视网膜病变视网膜OCT图像分割为n个区域，区域1、区域2、……区域n；

提取OCT图像的24个特征，具体方式包括：

步骤C1：设OCT图像中像素的坐标向量

从OCT图像标记区域中找到最上层上表面；从最上层上表面向下扫描，计算SF1下面像素到SF1的距离，作为一种距离特征；

OCT图像的横轴坐标x和纵轴坐标y作为另外两种特征，其中视网膜层的表面由z表示；

步骤C2：OCT图像

体素值也作为一种特征，使用各向异性滤波后的OCT图像

体素值也作为一种特征；

步骤C3：各向异性滤波后的OCT图像进行归一化操作，作为特征：

其中，

表示归一化后的图像，I_N,max归一化后的图像最大值，I_f,s表示归一化起始OCT图像值，I_f,r表示归一化亮度区间宽度；得到5个归一化图像特征；

步骤C4：计算图像

的每个像素三个特征值

其中σ_t表示尺度也即与图像

卷积的高斯函数的标准差，并计算亮层结构响应：

暗层结构响应：

其中，α、β表示对称性参数；基于公式(2)和公式(3)，按如下方式计算亮层和暗层结构的多尺度响应：

其中，σ_t,min,σ_t,max表示最小和最大尺度；示例随机森林训练与测试中，根据公式(2)-(4)计算OCT图像层结构响应图像特征14个；

提取特征后，使用分割标记区域和OCT图像特征训练随机森林分类器。

与一具体实施例中，在训练随机森林分类器时，将OCT图像的多个层结构分割为了8个分类标记区域，分区情况以及各层结构对应的上表面标号如下：区域1：玻璃体+脉络膜SF11’；区域2：神经纤维层SF1；区域3：神经节细胞层SF2；区域4：内丛状层SF3；区域5：内核层SF4；区域6：外丛状层SF5；区域7：外核层SF6+外界膜SF7+样区SF8；区域8：椭球区SF9+外光感受器节层10+交错区SF10’+视网膜色素上皮/布鲁赫SF11；

OCT图像层分割方法具体包括：

对待分割的OCT图像

使用高斯滤波器进行滤波，然后使用canny边缘检测算法找到SF1的初始形状，使用复合活性曲线算法找到最终的SF1；

找到最终的SF1后，提取OCT图像的24个特征，使用随机森林分类器，对OCT图像进行分类标记8个区域，从上向下找到标记区域的上表面作为视网膜8个区域的初始形状；

根据初始形状，使用复合活性曲线算法依次找到最终的区域2至区域7的上表面SF2-SF7；然后将区域7的上表面SF7向下平移10个像素作为区域8的上表面SF8的初始形状，使用复合活性曲线算法找到最终的SF8；同理，区域8的上表面SF8的最终形状向下平移若干像素点得到椭球区上表面SF9的初始形状，使用复合活性曲线算法找到最终的SF9；区域9的上表面SF9的最终形状向下平移若干像素点得到外光感受器节层上表面SF10的初始形状，使用复合活性曲线算法找到最终的SF10；

根据视网膜色素上皮/布鲁赫与脉络膜上表面SF11的初始形状，使用复合活性曲线算法找到最终的脉络膜的上表面SF11’；

使用阈值分割方法在SF10与SF11’之间分割出液体区域，即是最终的积液,计算积液的足迹；

向下平移5个像素作为SF10’，若在积液的足迹内，则找到积液的下表面更新SF10’；

SF11向上平移15个像素作为SF11，若在积液的足迹内，则找到积液的上表面更新SF11，且SF11上的点不能高于SF10’。

进一步地，确定区域1的上表面的初始形状具体包括：对待分割的OCT图像

使用高斯滤波器进行滤波，然后使用canny边缘检测算法找到区域1的上表面的的初始形状SF1。

进一步地，OCT图像分割为8个区域，区域1：玻璃体+脉络膜；区域2：神经纤维层；区域3：神经节细胞层；区域4：内丛状层；区域5：内核层；区域6：外丛状层；区域7：外核层+外界膜+样区；区域8：椭球区+外光感受器节层+交错区+视网膜色素上皮/布鲁赫。

进一步地，步骤C3中：I_f,s表示归一化起始OCT图像值，设为50、55、60、65、70；I_f,r表示归一化亮度区间宽度，设为120；

步骤C4中：α、β表示对称性参数，都设为0.25；

σ_t,min,σ_t,max表示最小和最大尺度，设为1和4,σ_t＝1,2,3,4)。

具体地，所述的预定的像素点为经验值。

与现有技术相比，本发明基于随机森林与复合活性曲线的OCT图像层分割方法具有以下优点：

首先，提取OCT图像的24个特征，基于该24个特征训练得到随机森林分类器，基于该随机森林分类器，得到OCT图像其中一些区域的上表面的初始形状，根据该一区域上表面初始形状使用复合活性算法得到该区域的上表面的最终形状。通过向上或向下平移预定像素点的方式得到剩余其他区域的上表面的初始形状，同样使用符合活性曲线算法得到该区域的上表面的最终形状。同理，当一个区域内包含多个层结构时，首先找到各层结构的初始形状，然后根据初始形状使用复合活性算法找到该层结构的上表面的最终形状。随机森林分类器训练精准，分割效果好精确，分割视网膜层与积液，克服现有的对于病变OCT图像分割算法识别率较低、分割效果较差等问题。

附图说明

图1为OCT的一个中心浆液性视网膜病变视网膜结构说明图像。

图2为OCT的一个的局部图像。

图3为复合活性曲线算法原理示意图，(a)OCT各向异性过滤图像与亮层亮层结构响应图像，黄色线为SF6的初始表面。(b)复合图G^H。(c)使用复合活性曲线算法分割得到精确表面SF6(红色线)叠加在OCT各向异性过滤图像与亮层亮层结构响应图像。

图4为表面分割对比图像，(a)原始OCT图像，(b)红色线为IF算法分割的表面而绿色线为专家手动分割的表面，(c)红色线为MGS算法分割的表面而绿色线为专家手动分割的表面，(d)红色线为随机森林方法找到的初始表面，(e)红色线为RFSGLW算法分割的表面而绿色线为专家手动分割的表面，(f)红色线为RFHLW算法分割的表面与σ_t,max＝4时的亮层结构响应图像，(g)为红色线为RFHLW算法分割的表面与σ_t,max＝2时的亮层结构响应图像，(h)为红色线为RFHLW算法分割的表面与而绿色线为专家手动分割的表面，(i)为RFHLW算法分割的表面三维显示。

图5为积液分割对比图像，(a)红色线为IF算法分割的积液而绿色线为专家手动分割的积液，(b)为GSGC算法分割的积液而绿色线为专家手动分割的积液，(c)为MGS算法分割的积液而绿色线为专家手动分割的积液，(d)为RFSGLW算法分割的积液而绿色线为专家手动分割的积液,(e)为RFHLW算法分割的积液而绿色线为专家手动分割的积液，(f)为RFHLW算法分割的积液三维显示；

说明：鉴于本发明的特殊性，为了更清楚地说明本发明的技术方案，在附图中部分采用了彩图的绘制。

具体实施方式

下面结合附图和实施例，对本发明的具体实施方式作进一步详细描述。以下实施例用于说明本发明，但不用来限制本发明的范围。

本发明基于随机森林与复合活性曲线的OCT图像层分割方法，包括：

训练随机森林分类器；在训练随机森林分类器时，将OCT图像的多个层结构分割为了n个分类标记区域，区域1、区域2、……区域n；提取了OCT图像中的m个特征；

确定区域1的上表面的初始形状，所述初始形状使用复合活性曲线算法得到区域1的上表面的最终形状SF1；

提取待分割OCT图像与训练随机森林分类器时相同的m个特征，使用随机森林分类器，对待分割OCT图像进行分类标记为n个区域，找到区域(n-i)、区域(n-j)……和/或区域(n-k)的各自上表面的初始形状，也即找到区域2至区域n中至少一个区域的上表面的初始形状。区域(n-i)、区域(n-j)……和/或区域(n-k)分别各自使用复合活性曲线算法得到区域(n-i)、区域(n-j)……和/或区域(n-k)的各自上表面SF(n-i)、SF(n-j)……和/或SF(n-k)的最终形状；

将所述的区域(n-i)、区域(n-j)……和/或区域(n-k)的各自上表面上表面SF(n-i)、SF(n-j)……和/或SF(n-k)的最终形状，中的某一个最终形状向上或向下移动若干像素点作为另一个剩余其他区域的上表面的初始形状，对该剩余其他区域的上表面的初始形状使用复合活性曲线算法得到该其他区域上表面的最终形状。

例如，将区域(n-i)上表面SF(n-i)的最终形状向上或向下平移a个像素点，作为区域(n-e)的初始形状，使用复合活性曲线算法得到区域(n-e)的上表面SFn(n-e)的最终形状。

同理，将区域(n-j)上表面SF(n-j)的最终形状向上或向下平移b个像素点，作为区域(n-f)的初始形状，使用复合活性曲线算法得到区域(n-f)的上表面SFn(n-f)的最终形状；……将区域(n-k)上表面SF(n-k)的最终形状向上或向下平移c个像素点，作为区域(n-g)的初始形状，使用复合活性曲线算法得到区域(n-g)的上表面SF(n-g)的最终形状。；

至于具体使用哪一个区域上表面的最终形状向上或向下移动若干像素点，根据实际情况作出具体地调整，移动多少个像素点，这个根据经验设定。i、j、k的具体取值，也根据具体经验设定。本发明中，i、j、k均为正整数。

还包括：判断一区域内是否包含OCT图像的至少两层结构，若包含至少两层结构，则首先得到其中一层结构的上表面的初始形状，使用随机森林活性算法得到该层结构的上表面的最终形状，然后该层结构的最终形状向上或向下平移预定像素点得到剩余其他一层结构的上表面的初始形状，基于该剩余其他层结构的上表面的初始形状，使用复合活性算法得到该剩余其他一层结构的上表面的最终形状。

其中，复合活性曲线算法具体包括：

步骤A1：对待分割的OCT图像使用各向异性滤波后的OCT图像

并计算亮层结构的多尺度响应响应

步骤A2：分别从

和

依次提取x-z平面图像，得到各自的二维图像I_2d和L_2d；提取出初始表面SF在当前切面中初始曲线

如果待分割曲面在OCT图像是由暗带到亮带交界面，则上下反转二维图像I_2d和L_2d；

如果待分割曲面在OCT图像是由亮带到暗带交界面，则不反转；

步骤A3：分别对二维图像I_2d和L_2d建立图G^F＝(V^F,E^F)和G^B＝(V^B,E^B)，其中V表示图的结点集合，E表示图的有向边集合；对于二维图像来说，图的结点是四个相邻像素交点，图的边是两个相邻像素p,q相接的边，有向边的右边像素表示内部、左边像素表示外部，按如下公式计算E^F和E^B向边e_ij的权值：

其中，f是特征函数，c_k是传递函数，ω_k是权值系数，n_l是边特征总数；n_k是传递函数总数；复合图G^H中的结点V^H与边E^H和图G^F中的结点V^F与边E^F、G^B中的结点V^B与边E^B一一对应，根据E^F和E^B计算复合图G^H＝(V^H,E^H)边的权值：

其中

表示点积；

将初始曲线

看成更新曲线

能量值E'←∞,E_lw←0，采用如下方式更新：

步骤B1：改变能量值E'←E_lw,E_lw←0，按间距lw_s从曲线

采样得到两个结点a_o和a_o+1，

a_x表示当前二维图像宽度；

步骤B2：在复合图G^H中，使用Dijkstra算法找到a_o到a_o+1的最小代价路径

其中，n₀表示当前路径途径结点总数；计算当前路径的代价值：

并更新代价值E_lw←E_lw+e_lw；代价路径

中等距采样得到lw_s结点更新曲线

步骤B3：更新采样步长lw_s←lw_s/2，如果|E_lw-E'|＞ΔE并且lw_s＞1，则返回步骤1；否则，平滑已更新曲线

如果输入图像上下反转，则把曲线

反转回去；最后，输出已更新曲线

本发明训练随机森林分类器的训练包括：

将将中心浆液性视网膜病变视网膜OCT图像分割为n个区域，区域1、区域2、……区域n；

提取OCT图像的24个特征，具体方式包括：

步骤C1：设OCT图像中像素的坐标向量

从OCT图像标记区域中找到最上层上表面；从最上层上表面向下扫描，计算SF1下面像素到SF1的距离，作为一种距离特征；

OCT图像的横轴坐标x和纵轴坐标y作为另外两种特征，其中视网膜层的表面由z表示；

步骤C2：OCT图像

体素值也作为一种特征，使用各向异性滤波后的OCT图像

体素值也作为一种特征；

步骤C3：各向异性滤波后的OCT图像进行归一化操作，作为特征：

其中，

表示归一化后的图像，I_N,max归一化后的图像最大值，I_f,s表示归一化起始OCT图像值，I_f,r表示归一化亮度区间宽度；得到5个归一化图像特征；

步骤C4：计算图像

的每个像素三个特征值

其中σ_t表示尺度也即与图像

卷积的高斯函数的标准差，并计算亮层结构响应：

暗层结构响应：

其中，α、β表示对称性参数；基于公式(2)和公式(3)，按如下方式计算亮层和暗层结构的多尺度响应：

其中，σ_t,min,σ_t,max表示最小和最大尺度；示例随机森林训练与测试中，根据公式(2)-(4)计算OCT图像层结构响应图像特征14个；

提取特征后，使用分割标记区域和OCT图像特征训练随机森林分类器。

下述各实施例中，以将OCT图像的多个层结构分割为了8个分类标记区域为例进行说明，分区情况以及各层结构对应的上表面标号如下：区域1：玻璃体+脉络膜SF11’；区域2：神经纤维层SF1；区域3：神经节细胞层SF2；区域4：内丛状层SF3；区域5：内核层SF4；区域6：外丛状层SF5；区域7：外核层SF6+外界膜SF7+样区SF8；区域8：椭球区SF9+外光感受器节层10+交错区SF10’+视网膜色素上皮/布鲁赫SF11。

在实际的应用过程中，还可以根据具体情况的不同将OCT图像分割为其他数量的标记区域，无论分区数量多少其各层上表面形状确定的原理是相同的，在此不再赘述。

实施例1

本实施例基于随机森林与复合活性曲线的OCT图像层分割方法，包括：

如图1所示，训练OCT图像特征训练随机森林分类器，在训练随机森林分类器时，将中心浆液性视网膜病变视网膜OCT图像分割为8个区域；区域1：神经纤维层；区域2：神经节细胞层；区域3：内丛状层；区域4：内核层；区域5：外丛状层；区域6：外核层+外界膜+样区；区域7：椭球区+外光感受器节层+交错区+视网膜色素上皮/布鲁赫；区域8(class 8)：玻璃体+脉络膜；区域1的上表面层Surface1、区域2的上表面Surface2……区域8的上表面，Surface1简称SF1，同理Surface2简称SF2……Surface8简称SF8，具体参见图1所示。

OCT图像层分割具体方法，包括：

对待分割的OCT图像

使用高斯滤波器进行滤波，然后使用canny边缘检测算法找到SF1的初始形状。然后，使用复合活性曲线算法找到精确的SF1。

提取OCT图像的24个特征，使用随机森林分类器，对OCT图像进行分类标记8个区域，从上向下找到标记区域的上表面作为视网膜8个区域的初始形状；

根据初始形状，使用复合活性曲线算法依次找到最终的SF2至SF7；然后将SF7向下平移10个像素作为SF8的初始形状，使用复合活性曲线算法找到最终的SF8；类似地，找到最终的SF9、SF10；

视网膜色素上皮/布鲁赫与脉络膜的初始形状使用复合活性曲线算法找到最终的SF11’；

使用阈值分割方法在SF10与SF11’之间分割出液体区域，即是最终的积液，计算积液的足迹；

SF10向下平移5个像素作为SF10’，若在积液的足迹内，则找到积液的下表面更新SF10’；

SF11向上平移15个像素作为SF11，若在积液的足迹内，则找到积液的上表面更新SF11，且SF11上的点不能高于SF10’。

实施例2

本实施例基于随机森林与复合活性曲线的OCT图像层分割方法，在实施例1的基础上，得到OCT图像特征训练随机森林分类器具体包括：

将中心浆液性视网膜病变视网膜OCT图像分割为8个区域；区域1：神经纤维层；区域2：神经节细胞层；区域3：内丛状层；区域4：内核层；区域5：外丛状层；区域6：外核层+外界膜+样区；区域7：椭球区+外光感受器节层+交错区+视网膜色素上皮/布鲁赫；区域8(class8)：玻璃体+脉络膜；区域1的上表面层Surface1、区域2的上表面Surface2……区域8的上表面，Surface1简称SF1，同理Surface2简称SF2……Surface8简称SF8，具体参见图1所示该具体划分不限于实施例1和本实施例2中所述的分类规则，可以根据具体需要进行调整。

提取OCT图像的24个特征，具体方式如下：

步骤C1：设OCT图像中像素的坐标向量

从OCT图像标记区域中找到最上层上表面；从最上层上表面向下扫描，计算SF1下面像素到SF1的距离，作为一种距离特征；OCT图像的横轴坐标x和纵轴坐标y作为另外两种特征，其中视网膜层的表面SF可以由z表示；图2说明了OCT图像坐标系、SF6和OCT图像三个平面切片图；

步骤C2：OCT图像

也作为一种特征，使用各向异性滤波后的OCT图像

也作为一种特征；

步骤C3：各向异性滤波后的OCT图像进行归一化操作，也可以作为特征：

其中，

表示归一化后的图像，I_N,max归一化后的图像最大值(示例中设为255)，I_f,s表示归一化起始OCT图像值(示例中设为50、55、60、65、70)，I_f,r表示归一化亮度区间宽度(示例中设为120)；示例中得到5个归一化图像特征；

步骤C4：计算图像

的每个像素三个特征值

其中σ_t表示尺度(与图像

卷积的高斯函数的标准差)，并计算亮层结构响应：

暗层结构响应：

其中，α、β表示对称性参数(示例中都设为0.25)；基于公式(2)和公式(3)，可按如下方式计算亮层和暗层结构的多尺度响应：

其中，σ_t,min,σ_t,max表示最小和最大尺度；示例随机森林训练与测试中，根据公式(2)-(4)计算OCT图像层结构响应图像特征14个；

提取特征后，基于OCT图像分割标记和OCT图像特征训练随机森林分类器。

本发明所采用的OCT图像数据来自江苏省人民医院，其中6个数据作为训练数据，48个数据作为分析此发明方法分割结果准确性的比较。通过真阳性分数TPF、假阳性分数FPF、DSC系数(Dice similarity coefficient)来直观统计积液分割结果与金标准的重合度，即分割的准确性，定义为：

|·|表示体素个数，F_r表示专家手动标记为积液，B_r表示专家手动标记为非积液，F_a表示发明方法分割结果为积液。将此发明方法(RFHLW)与现有方法Iowa参考算法(IF)、图搜素图割算法(GSGC)、多分辨率图搜素算法(MGS)以及随机森林+活性曲线算法相比较积液分割准确性，如表1所示，做配对t-test测试，如表2所示，p值<0.05表示两种方法显著差异性大。

表1

表2

通过计算专家手动标记表面与此发明方法分割表面之间的无符号距离和符号距离，评估此发明方法表面分割的精度。将此发明方法与现有方法Iowa参考算法(IF)、多分辨率图搜素算法(MGS)以及随机森林+活性曲线算法(RFSGLW)相比较，做配对t-test测试，p值<0.05表示两种方法显著差异性大。无符号距离对比如表3所示，符号距离如表5所示，表4为此发明方法无符号距离与其余方法的t-test测试p值，表6为此发明方法符号距离与其余方法的t-test测试p值.

表3

表4

表5

表6

以上所述仅是本发明的优选实施方式，并不用于限制本发明，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明技术原理的前提下，还可以做出若干改进和变型，这些改进和变型也应视为本发明的保护范围。

Claims (7)

1.一种基于随机森林与复合活性曲线的OCT图像层分割方法，其特征在于，包括：

训练随机森林分类器；在训练随机森林分类器时，将OCT图像的多个层结构分割为了n个分类标记区域，各区域分别具有上表面，各区域的上表面层由上至下依次为SF1、SF2、SF3……SFn；提取了OCT图像中的m个特征；

确定待分割OCT图像的上表面SF1的初始形状，使用复合活性曲线算法得到上表面SF1的最终形状；

提取待分割OCT图像与训练随机森林分类器时相同的m个特征，使用随机森林分类器，对待分割OCT图像进行分类标记为n个区域，找到区域2至区域n中N个区域上表面的初始形状，基于各区域上表面的初始形状分别各自使用复合活性曲线算法，得到各区域上表面的最终形状，其中n、N均为正整数，N≤n-1；

当N＜n-1时，还包括：将一区域的上表面的最终形状向上或向下平移预定像素点得到剩余其他一区域的上表面的初始形状，基于该剩余其他一区域的上表面的初始形状，使用复合活性算法得到该剩余其他一区域的上表面的最终形状；

还包括：判断一区域内是否包含OCT图像的至少两层结构，若包含至少两层结构，则首先得到其中一层结构的上表面的初始形状，使用随机森林活性算法得到该层结构的上表面的最终形状，然后该层结构的最终形状向上或向下平移预定像素点得到剩余其他一层结构的上表面的初始形状，基于该剩余其他一层结构的上表面的初始形状，使用复合活性算法得到该剩余其他一层结构的上表面的最终形状；

其中，复合活性曲线算法具体包括：

步骤A1：对待分割的OCT图像使用各向异性滤波后的OCT图像

并计算亮层结构的多尺度响应响应

其中，

表示体素的坐标，σ_t表示尺度；

步骤A2：分别从

和

依次提取x-z平面图像，得到各自的二维图像I_2d和L_2d；提取出初始表面SF在当前切面中初始曲线

如果待分割曲面在OCT图像是由暗带到亮带交界面，则上下反转二维图像I_2d和L_2d；

如果待分割曲面在OCT图像是由亮带到暗带交界面，则不反转；

步骤A3：分别对二维图像I_2d和L_2d建立图G^F＝(V^F，E^F)和G^B＝(V^B，E^B)，

其中，V表示图的结点集合，E表示图的有向边集合；对于二维图像来说，图的结点是四个相邻像素交点，图的边是两个相邻像素p，q相接的边，有向边的右边像素表示内部、左边像素表示外部，

按如下公式计算E^F和E^B向边e_ij的权值：

其中，f是特征函数，c_k是传递函数，ω_k是权值系数，n_l是边特征总数；n_k是传递函数总数；复合图G^H中的结点V^H与边E^H和图G^F中的结点V^F与边E^F、G^B中的结点V^B与边E^B一一对应，根据E^F和E^B计算复合图G^H＝(V^H,E^H)边的权值：

其中

表示点积；

将初始曲线

看成更新曲线

能量值E'←∞,E_lw←0，采用如下方式更新：

步骤B1：改变能量值E'←E_lw,E_lw←0，按间距lw_s从曲线

采样得到两个结点a_o和a_o+1，

a_x表示当前二维图像宽度；

步骤B2：在复合图G^H中，使用Dijkstra算法找到a_o到a_o+1的最小代价路径

其中，n₀表示当前路径途径结点总数；计算当前路径的代价值：

并更新代价值E_lw←E_lw+e_lw；代价路径

中等距采样得到lw_s结点更新曲线

步骤B3：更新采样步长lw_s←lw_s/2，如果|E_lw-E'|＞ΔE并且lw_s＞1，则返回步骤B1；否则，平滑已更新曲线

如果输入图像上下反转，则把曲线

反转回去；最后，输出已更新曲线

2.根据权利要求1所述的基于随机森林与复合活性曲线的OCT图像层分割方法，其特征在于，训练随机森林分类器的训练包括：

将将中心浆液性视网膜病变视网膜OCT图像分割为n个区域，区域1、区域2、……区域n；

提取OCT图像的24个特征，具体方式包括：

步骤C1：设OCT图像中体素的坐标向量

从OCT图像标记区域中找到最上层上表面；从最上层上表面向下扫描，计算SF1下面像素到SF1 的距离，作为一种距离特征；

OCT图像的横轴坐标x和纵轴坐标y作为另外两种特征，其中视网膜层的表面由z表示；

步骤C2：OCT图像

的体素值也作为一种特征，使用各向异性滤波后的OCT图像

的体素值也作为一种特征；

步骤C3：各向异性滤波后的OCT图像进行归一化操作，作为特征：

其中，

表示归一化后的图像，I_N,max归一化后的图像最大值，I_f,s表示归一化起始OCT图像值，I_f,r表示归一化亮度区间宽度；得到5个归一化图像特征；

步骤C4：计算图像

的每个体素三个特征值

其中σ_t表示尺度也即高斯函数的标准差，并计算亮层结构响应：

暗层结构响应：

其中，α、β表示对称性参数；基于公式

(2)和公式(3)，按如下方式计算亮层和暗层结构多尺度响应：

其中，σ_t，min,σ_t,max表示最小和最大尺度；示例随机森林训练与测试中，根据公式(2)-(4)计算OCT图像层结构响应图像特征14个；

提取特征后，使用分割标记区域和OCT图像特征训练随机森林分类器。

3.根据权利要求2所述的基于随机森林与复合活性曲线的OCT图像层分割方法，其特征在于，

在训练随机森林分类器时，将OCT图像的多个层结构分割为了8个分类标记区域，分区情况以及各层结构对应的上表面标号如下：区域1：神经纤维层SF1；区域2：神经节细胞层SF2；区域3：内丛状层SF3；区域4：内核层SF4；区域5：外丛状层SF5；区域6：外核层SF6+外界膜SF7+样区SF8；区域7：椭球区SF9+外光感受器节层SF10+交错区SF10’+视网膜色素上皮/布鲁赫SF11；区域8：玻璃体+脉络膜SF11’；

OCT图像层分割方法具体包括：

对待分割的OCT图像

使用高斯滤波器进行滤波，然后使用canny边缘检测算法找到SF1的初始形状，使用复合活性曲线算法找到最终的SF1；

找到最终的SF1后，提取OCT图像的24个特征，使用随机森林分类器，对OCT图像进行分类标记8个区域，从上向下找到标记区域的上表面作为视网膜8个区域的初始形状；

根据初始形状，使用复合活性曲线算法依次找到最终的区域2至区域7的上表面SF2-SF7；然后将区域7的上表面SF7向下平移10 个像素作为区域8的上表面SF8的初始形状，使用复合活性曲线算法找到最终的SF8；同理，区域8的上表面SF8的最终形状向下平移若干像素点得到椭球区上表面SF9的初始形状，使用复合活性曲线算法找到最终的SF9；区域9的上表面SF9的最终形状向下平移若干像素点得到外光感受器节层上表面SF10的初始形状，使用复合活性曲线算法找到最终的SF10；

根据视网膜色素上皮/布鲁赫与脉络膜上表面SF11的初始形状，使用复合活性曲线算法找到最终的脉络膜的上表面SF11’；

使用阈值分割方法在SF10与SF11’之间分割出液体区域，即是最终的积液,计算积液的足迹；

向下平移5个像素作为SF10’，若在积液的足迹内，则找到积液的下表面更新SF10’；

SF11向上平移15个像素作为SF11，若在积液的足迹内，则找到积液的上表面更新SF11，且SF11上的点不能高于SF10’。

4.根据权利要求1所述的基于随机森林与复合活性曲线的OCT图像层分割方法，其特征在于，确定区域1的上表面SF1的初始形状具体包括：对待分割的OCT图像

使用高斯滤波器进行滤波，然后使用canny边缘检测算法找到区域1的上表面的初始形状SF1。

5.根据权利要求1所述的基于随机森林与复合活性曲线的OCT图像层分割方法，其特征在于，OCT图像分割为8个区域，区域1：玻璃体+脉络膜；区域2：神经纤维层；区域3：神经节细胞层；区域4：内丛状层；区域5：内核层；区域6：外丛状层；区域7：外核层+外界膜+样区；区域8：椭球区+外光感受器节层+交错区+视网膜色素上皮/布鲁赫。

6.根据权利要求2所述的基于随机森林与复合活性曲线的OCT图像层分割方法，其特征在于，步骤C3中：I_f,s表示归一化起始OCT图像值，设为50、55、60、65、70；I_f,r表示归一化亮度区间宽度，设为120；

步骤C4中：α、β表示对称性参数，都设为0.25；

σ_t,min,σ_t,max表示最小和最大尺度，设为1和4,σ_t＝1，4。

7.根据权利要求1所述的基于随机森林与复合活性曲线的OCT图像层分割方法，其特征在于，所述的预定的像素点为经验值。

Images