陆圣陶+陈强+牛四杰+王玉萍

摘 要： 频域光学断层相干扫描（SDOCT）图像针对视网膜中的RNFL层的精确检测，对于一些眼科疾病诊断有很好的参考价值，需要精确分割该层的上表面ILM边界。给出一种3D图搜索与梯度弱化方式相结合的分割ILM表面的方法，提出一种梯度弱化的方法，弱化了图像中其他高对比度区域的梯度，降低其他层对ILM表面分割的干扰，从而在图像中仅有在ILM表面附近较大梯度的前提下，利用经典三维图搜索方法实现对ILM表面的精确分割。分割过程中采用多尺度的分割方法，降低了三维图搜索算法本身时间复杂度和空间复杂度。采用斯坦福大学的25对眼睛的频域光学相干断层（SDOCT）图像，共50组3D图像体数据实验，结合ITK?SNAP软件的半自动分割结果进行实验对比。实验结果表明给出的方法对50组数据中的43组数据能够实现ILM表面的快速精确分割。

关键词： 3D图搜索； 梯度弱化； 多尺度分割； ILM表面分割

中图分类号： TN919?34； TP391 文献标识码： A 文章编号： 1004?373X（2015）04?0106?05

0 引 言

OCT图像层分割在常见的眼科疾病中（青光眼、AMD（Age?Related Macular Degeneration，老年黄斑性变性））对于病变区域的诊断有很高的参考价值 [1]。图1给出了一帧黄斑中央凹附近区域单帧的图像分层示意图。

近年来针对视网膜分层的方法主要包含如下几类。文献[2]提出基于区域增长的分割方法，该方法不能保证所求表面全局最优；Kajic等人给出了一种基于图像纹理和形状统计信息的边缘特征提取分割算法，医学图像的纹理特征提取相对比较困难[3]； 文献[4]采用结合形态学的边缘提取方法，该方法针对部分边缘不明显的医学图像效果不甚好，没有利用帧间信息；文献[5]给出了基于水平集方法分割OCT图像，其对应的符号距离函数的计算量较大。以上几类分割方法均以OCT图像的一个B?Scan（一帧切片图像）为单位进行图像分割，没有充分考虑连续图像帧与帧之间相对连续关系；目前基于图论的分割方法可分成2D和3D两类。2D的图论算法的图像分割单位针对单帧图像。文献[6]采用图搜索的方法对OCT视网膜图像进行层分割，利用Dijkstra算法搜索出最短路径，逐步区域限制分割出各个视网膜层，其针对一些边缘不明显图像也需要做部分分割后处理操作。基于3D图论的方法在近年来得到广泛的应用。3D图搜索采用构图的方法，给出了多表面同时分割的理念并引入3D子图间限制，最终转换成求解一个4D结构图的最优解问题，目前Dufour等人在3D图搜索基础上又给出了先验信息模型的概念，并结合先验信息模型给出了软限制构图的思想，但依然需要层与层之间约束构图同时分割求解，不能做到简单，单独分割某层[7]。如图1所示。其中表面标记从上到下分别为内界膜（Internal Limiting Membrane，ILM），INL?OPL（Inner nuclear layer? Outer plexiform layer），OPL?ONL（Outer plexiform layer? Outer nuclear layer），IS?OS（photoreceptor innerand outer segments），OS?RPE（outer segments?retinal pigment epithelium）和BM（Bruch's membrane）。

直接使用3D图搜索算法求解较慢，并且难以单独分割某一层。本文在结合3D图搜索方法基础上，引入一种梯度弱化的方法，并将该方法应用到ILM表面的分割问题上。该方法优势是，第一，能够单独分割ILM层。直接单独分割ILM层可能会产生误分割，常采用3层同时分割才可分出ILM层，当仅需分割ILM层时不甚合理；第二，利用梯度弱化后结合3D图搜索分割ILM层会加快图搜索的速度，原因是仅保持ILM表面附近较大梯度，其他区域梯度较小且相对一致，顶点权重w的计算方式导致构图过程很多顶点权重为0，在将3D图搜索转化成求解最小割问题后，这些点会排除在s?t图以外，大大减少s?t图中的边数，加快求解最小割速度。

1 方 法

本节给出梯度弱化方法的概念及算法实现。算法实现流程如图2所示。

1.1 3D图搜索算法模型

3D图搜索算法能量函数分成权值能量部分项和平滑约束能量部分项。构造图中的边分成两种类型有向边：列内偏序关系边和列间限制关系边。假设每列当中上下相邻的两个顶点坐标形式为V（x，y，z），V（x，y，z-1），则列内偏序关系边可表示为如下形式：

[Edgeintra={V（x，y，z），V（x，y，z-1）|z≥1}] （1）

对于列间限制关系边，假设在x和y方向的高度差限制分别为Δx和Δy，相邻列的点构造如下列间限制边：

[Edgeinter = {V（x，y，z），V（x+1，y，max（z-Δx，0）） V（x，y，z），V（x-1，y，max（z-Δx，0）） V（x，y，z），V（x，y-1，max（z-Δy，0）） V（x，y，z），V（x，y+1，max（z-Δy，0））}] （2）

这两种关系边的权值设定为∞，其实际意义在随后构造s?t图时该条路径的流容量无穷大。列间限制为简单起见采用固定的Δx和Δy，实验过程令Δx=Δy=3。这两种限制保证了帧与帧图像间的限制关系。最终求最优曲面可转化成求图的最小闭集。可证明该最小闭集的问题可转化成求图的s?t割问题。具体基础理论内容可参考相应文献。

1.2 梯度弱化算法

对图像进行双边滤波的医学图像去噪[7]，并进行灰度填充预处理去噪之后采用梯度弱化算法。梯度弱化是指将所感兴趣之外区域的梯度抑制，降低其他部分对所分割区域的干扰，具体步骤如下。

1.2.1 背景估计

针对每一帧OCT图像背景和目标阈值分割的阈值通过式（3）确定：

[ti=max0≤k≤255h（Iik）+Δg， i=1，2，…，128] （3）

式中：t表示需要获取的阈值；i表示体数据的帧数；k表示灰度级；函数[h（k）]计算灰度级k的像素个数；Δg表示容许的灰度差异，可取一个较小的常量，比如Δg=10。通过获取得到阈值之后对图像进行阈值分割即可先把背景区域分离出来。

1.2.2 视网膜区域估计

通过视网膜区域估计得到感兴趣区域。如上小节所述，得到的目标区域不一定是视网膜区域，也有可能是其他部分的较亮的条纹或者亮斑带来的干扰。通过取最大连通区域来消除干扰，计算步骤如下：

（1） 对二值图像连通区域标记。假设二值图像连通区域个数为n，对应的连通区域分别为[p1，p2，…，pn]，[j∈pipij]统计连通区域p内目标像素个数。二值图像中两个目标像素属于同一个4?连通区域条件是在4?邻域条件下可达。对每帧二值图像可通过式（7）计算最大连通分量：

[pmax=maxi=1，2，...，n（j∈pipij）] （4）

（2） 在（1）基础上，获取每一列最大连通区域，用来消除视网膜上方玻璃体条带状物质带来的干扰。二值图像仅保留最大连通分量[pmax]。假设二值图像第i列连通分量个数为ni，对应的连通区域分别为[ci1 ， ci2，… ，cini]，并且[j∈cimcimj]表示统计连通区域[cim]内目标像素个数。通过式（5）计算每一列的最大连通分量（col表示图像的列数）：

[cimax=maxm=1，2，...，ni（j∈cimcimj） ， i=1，2，…，col] （5）

在经过以上步骤后得到视网膜目标区域在OCT图像中的位置。

1.2.3 目标区域腐蚀及梯度弱化

目标区域腐蚀目的是获得非感兴趣区域，并将该区域梯度弱化，主要思想是对得到的视网膜估计的感兴趣区域进行形态学腐蚀得到非感兴趣区域。假设目标区域的单帧灰度图像为I，其对应的区域估计后的二值图像为Ibw。首先通过式（6）对二值图像进行腐蚀，SE表示腐蚀操作的结构元素，这里采用半径为3的圆作为结构元素，*为卷积操作。腐蚀操作后的目标区域会收缩变小，从而能够保证腐蚀后保留出ILM表面附近区域：

[Mdata=Ibw? SE] （6）

[Iout=I*-1-2-1000121?（Mdata）] （7）

将腐蚀之后的目标区域内的梯度通过式（7）去除，其中[Mdata]表示对图像进行取反操作。去除后保证在输出图像中，仅有ILM表面附近的强梯度区域得以保留，最终得到的较大梯度的感兴趣区域仅为ILM附近区域。

1.3 ILM表面分割

以弱化后的梯度图像来计算每个体素的权值，并对整张图利用3D的图搜索算法进行操作。采用文献[9]使用的最大流最小割算法计算最小割。具体步骤如下：

（1） 对图像4倍下采样，利用3D图搜索方法进行分割求解，得到下采样ILM层分割结果；

（2） 利用对下采样ILM层分割结果采用文献[10]使用的薄板样条函数曲线拟合得到粗分割结果；

（3） 对原始图像ILM层细分割，在限制范围内采用3D图搜索算法，最终得到细分割结果。

2 实验结果和分析

本文采用了斯坦福大学25对眼睛的SD?OCT图像进行实验，图像体数据分辨率为512×1 024×128。每只眼睛分成OD（右），OS（左）两组图像，共50组图像，利用Matlab脚本编程与C++混编完成算法过程实现。使用ITK?SNAP软件的snake算法进行半自动分割ILM表面，并利用该软件本身自带手动分割功能进行调整后作为标准，与本文算法得出的结果进行对比，同时也对50组图像直接采用3D图搜索算法分割ILM表面，给出了本文算法和直接利用3D图搜索分割ILM层的算法的实验结果精度对比。

2.1 定性分析

实验对不同的SD?OCT图像进行分割，并在图3中给出了单帧正常图像的分割结果。在50组数据集中主要包括两种特殊情况，一是在正常图像里中央凹区域不明显，二是分割病变图像，在ILM层附近主要可能的病变为下积液。图3给出这两种特别情况的分割结果。图3（c） 比图3（b）分割更加精确，对图像灰度填充降低了下积液对ILM分割的干扰，不会分割到下积液里去。图3（f）对中央凹区域梯度不明显分割效果较精确，原因是结合3D图搜索本身的列间限制和薄板样条插值保证分割结果光滑；直接3D图搜索加薄板样条插值可能由于中央凹区域较暗，该情况会有可能发生误分割，图3（e）给出发生误分割的情况。

2.2 定量分析

表1给出误差的平均结果对比（误差的平均结果不包括误分割的数据在内），加粗数据为误分割。实验采用ITK软件进行半自动分割得出的ILM表面为基准，分析了本文算法得出ILM表面的位置与基准位置的误差，并且在图4中分别给出了OD（右）和OS（左）的对比数据统计结果（由于误分割均值和标准差较大，为便于观察对比结果，将误分割的均值换成了一个相对较小的值8，方差设为0，如图4所示）。

结合图4和表1数据，不包含误分割情况，直接采用3D图搜索误差平均结果为OD（右）：（2.68±2.50）μm，OS（左）：（2.94±2.98）μm，本文算法25组OD（右）误差的平均结果为（2.55±2.35） μm，OS（左）为（2.71±2.73） μm，比直接采用3D图搜索精度要高，而且避免了误分割情况的发生。

2.3 分析与讨论

（1） 分割精度。由表1可以看出，50组数据中，有47组分割结果比直接利用3D图搜索算法分割ILM层精度要高，原因是在获得梯度弱化后的梯度图像前，算法进一步去除了分割过程中产生的部分干扰，能够在一定程度上提高分割精度。须注意加粗的异常分割结果，产生原因是直接利用3D图搜索算法分割ILM层，由于其他层可能梯度较大，而求解过程为获得全局最优解，在不改变权重的情况下会造成误分割，如图5所示。本文算法由于弱化了其他层，避免了误分割的发生。

（2） 时间复杂度。直接采用3D图搜索算法分割128帧图像需要300 s以上的时间；而本文算法所有数据都能在70 s以内分割128帧图像，平均单张分割不到1 s。原因是由于梯度弱化后在构造s?t图时由于很多权重为0的顶点不会与源点s和汇点t有边，从而降低了求解图割的构图规模，从而降低了计算的时间复杂度。

（3） 算法鲁棒性。能够精确分割ILM边界，并对可能出现的中央凹不明显以及存在下积液孔洞的特殊情况鲁棒性较好，如图5所示。

3 结 语

本文引入梯度弱化的概念分割SDOCT图像中的ILM层，利用3D图搜索算法考虑到图像帧间的互信息的优势，并通过梯度弱化过程分割ILM表面，降低了其他区域对该表面分割的干扰，提高了分割的精度和图搜索的速度。但对于图像体素的权值函数应该结合其他更多的信息。另外，对于3D图搜索仅采用比较基本的思想，表面限制采用单一的硬限制，今后可考虑将3D图搜索中的多层同时分割引入梯度弱化思想，考虑构造比较完善的权值函数。

参考文献

[1] 孙延奎.光学相干层析医学图像处理及其应用[J].光学精密工程，2014（4）：1086?1104.

[2] 徐奇，常英，李文彬，等.基于区域生长的OCT图像分层算法[J].信息技术，2013（4）：136?140.

[3] KAJIC V， ESMAEELPOUR M， POVAZAY B， et al. Automated choroidal segmentation of 1060 nm OCT in healthy and pathologic eyes using a statistical mode [J]. Biomed Opt Express， 2012， 3（1）： 86?103.

[4] 张博闻，田小林，孙延奎.基于改进的数学形态学的OCT图像快速边缘提取算法[J].计算机科学，2013（z1）：173?175.

[5] 杨平，彭清，刘维平，等.一种眼底黄斑水肿OCT图像分割方法[J].生物医学工程学杂志，2011（5）：1001?1006.

[6] 胡超.基于OCT图像玻璃疣的自动检测与分割[D].南京：南京理工大学，2013.

[7] DUFOUR P A， CEKLIC L， ABDILLAHI H. Graph?based multi?surface segmentation of OCT data using trained hard and soft constraints [J]. IEEE Transactions on Medical Imaging， 2013， 32（3）： 531?542.

[8] 张聚，王陈，程芸.小波与双边滤波的医学超声图像去噪[J].中国图象图形学报，2014（1）：126?132.

[9] 张铃.动态网络上最大流概念及其性质的研究[J].模式识别与人工智能，2013（7）：609?614.

[10] 任伟建，张宏鑫.基于薄板样条函数的地图矢量化处理[J].系统仿真技术，2012（3）：187?191.