CN108961261B - 一种基于空间连续性约束的视盘区域oct图像层次分割方法 - Google Patents

Abstract

本发明公开了一种基于空间连续性约束的视盘区域OCT图像层次分割方法，该方法首先使用BM3D算法去除OCT图像中的散斑噪声；之后使用基于模糊C均值和主动轮廓的ROI分割方法将图像中多个高反射区相互分离；然后使用ROI区域图像对血管阴影进行定位；使用A‑scan分割算法依次对所述ROI区域进行分割，利用血管阴影区域对每幅图像的初步分割结果进行修正；最后使用空间连续性约束优化修正后的分割结果，获得ILM、IS‑OS、BM的分割边界；该方法是一种行之有效的视盘区域OCT图像的层次分割方法，分割准确度较高，并对散斑噪声、血管阴影具有一定的鲁棒性。

Description

一种基于空间连续性约束的视盘区域OCT图像层次分割方法

技术领域

本发明属于医学图像处理和图像分割技术领域，特别涉及一种基于空间连续性约束的视盘区域OCT图像层次分割方法。

背景技术

频域光学相干断层扫描(Spectial-Domain Optical Coherence Tomography,SD-OCT)成像技术自引入以来，由于其成像过程中对拍摄对象具有非侵害性的特点，被广泛用于现代眼科医学影像诊断中。相比较于眼底照相机、超声波等传统医学影像技术而言，SD-OCT能够获得更加丰富的生理组织的三维结构信息，为医生的临床诊断提供更加全面化的辅助。

常见的眼科疾病，如青光眼、老年性黄斑变性(Age-related MacularDegeneration,AMD)等，大多属于视网膜疾病。这些眼科疾病发生时，大多数患者的视网膜生理结构层次会发生异常变化，特别是视网膜ONH(Optical Nerve Head)区域和黄斑区域。通过对患者病例进行分析发现，视网膜中常见的异常变化主要有：视网膜RNFL(RetinalNerve Fiber Layer)层厚度变薄、ONH区域视杯视盘面积变小、杯盘比(Cup/Disc Ratio，CDR)变大等。因此，如何对这些视网膜结构参数进行准确测量，是眼科疾病计算机辅助诊断的关键步骤。同时，这些结构参数的测量涉及到对视网膜OCT图像的层次分割。由于目前的OCT系统能够获得视网膜区域连续的断层图像，导致OCT图像数据量巨大。如果采用人工的方式对这些图像进行分割，不仅费时费力而且无法保证多张图像分割结果的一致性，导致参数测量的不准确。因此，如何借助计算机辅助技术快速有效地计算出视网膜生理层厚度以及对视网膜结构参数进行测量，是亟需解决的一个问题。针对这个问题，国内外许多研究人员对OCT图像层次分割算法进行了研究。目前所公开的OCT层次分割算法可以归为五个不同的类别，分别是基于A-scan、B-scan、Active contours、Artificial intelligence和3Dgraphs等。

(1)基于A-scan的分割算法：Hee首次提出基于A-scan的方法对视网膜RNFL的厚度进行测量。该方法利用图像中每个A-scan上的像素值的变化进行一维边缘检测，通过检测出的边缘点来确定RNFL层的上下界面，之后测量RNFL层的厚度。在A-scan分割的基础上，Koozekanani结合视网膜本身的结构特性，提出了一些有关分割边界的假设，例如，正常视网膜的边界是平滑且连续的，每个A-scan与每条边界仅有一个交点。该方法在这些假设的基础上，使用马尔科夫序列(Markov sequence)对相邻的分割结果进行建模，最终获得较为准确的分割边界。Shahid和Ishikawa提出了基于图像像素值的变化分割黄斑区域视网膜层次的算法，该方法中，Ishikawa提出了黄斑中心凹处的视网膜厚度能够作为眼科疾病的诊断依据；同时，该方法在预处理过程中采用了A-scan之间的配准和直方图均衡化的操作，减少了噪声对分割结果的影响。经过分析发现，这些基于A-scan的分割算法大多没有结合相邻A-scan之间的相似性以及相邻图像之间的位置对应关系，其分割结果中会出现生理层错分的情况，把目标生理层的一部分错误分割到上下临近的生理层上，导致最终的分割准确度较差。同时这些方法不能够很好地处理含有大量斑点噪声的图像，导致分割准确率较差。

(2)基于B-scan的分割算法：Boyer等在Koozekanani方法基础上改进，使用抛物线模型和马尔科夫模型对分割结果进行优化，提取了一些视网膜的结构性参数用作眼科疾病临床诊断的依据；其中，CDR被阐述可以作为眼科疾病诊断的依据。之后，Baroni等通过在分割过程中最大化边界似然函数来获得准确的分割结果。该方法在分割步骤中，首先根据边界附近像素的梯度特性构建边界似然函数，之后在A-scan分割结果中选择能够最大化该似然函数的分割点作为最终分割结果中的边界点，最终所有符合条件的点组成目标边界。Tan等利用动态规划方法和图像的二维梯度信息较为准确地提取到视网膜层次的厚度信息，同时展示了视网膜中不同层次厚度之间的关系，例如RNFL、神经节细胞层(Ganglion CellLayer,GCL)和内丛状层(Inner Plexiform Layer,IPL)。最近，Kajic等在预处理步骤中使用基于双树复小波变换(Dual-tree complex wavelet,DTCW)降噪算法对OCT图像进行处理，能够较好地去除图像中的斑点噪声；之后使用基于图像纹理和视网膜形状的统计模型对OCT图像进行分割，能够准确分割得到视网膜的八个层次。相比较于A-scan的方法而言，B-scan方法能够在预处理步骤中更好地去除图像中的斑点噪声，同时在分割过程中结合邻域图像的相似性，可以获得更加准确的分割结果。

(3)基于Active contours的分割算法：Active contours模型通过最小化当前轮廓的内部区域和外部区域的能量和来对图像中的边界进行分割。Zhu等针对ONH区域的图像，提出了基于主动轮廓的Floatingcanvas方法用于视网膜内三维层次的分割，该方法可以获得相对平滑的边界面，但在边界附近梯度较低的情况下会出现分割错误。Yazdanpanah等结合有关视网膜层次形状的先验信息，利用Active contours方法对图像进行分割，该方法能够有效地克服血管阴影造成的影响，但在视网膜形状不规则的情况下分割效果较差。Ghorbel等提出了一种结合主动轮廓和马尔科夫随机场理论的分割方法，该方法能够准确地将视网膜分为八个层次，但图像中的血管阴影对分割结果的准确率影响较大。这类方法在图像对比度高、视网膜层次边界平滑的图像上分割准确度高，然而在有血管阴影和视网膜形状不规则的情况下效果较差。

(4)基于Artificial intelligence的分割算法：这类方法的基本原理是训练一个二分类器，之后使用这个分类器对像素点进行分类，判断其是否属于视网膜层次边界，该分类器的训练过程需要以专家手工分割的数据作为金标准。Fuller等人设计了一种使用多分辨率图像的支持向量机(Support Vector Machines，SVMs)进行视网膜OCT图像的层次分割与厚度测量的方法；相比较于其它分割方法而言，该方法的分割准确度略低且有8％的厚度测量误差。之后，Mayer等人使用模糊C均值聚类算法(Fuzzy C-means clustering，FCM)对来自5个正常人和7个青光眼患者的视网膜图像进行RNFL厚度的测量，该方法的RNFL层分割结果相对于金标准而言，分别有97％上边界点和74％下边界点落在与金标准相差两个像素的范围内。王玉萍等通过使用支持向量机、极限学习机等对视网膜RPE层进行分割，并通过与杯盘比的结合探讨了青光眼的诊断标准。这类方法的分割结果准确度依赖于训练集的选择，对数据要求比较高。

(5)基于3D graph的分割算法：基于图论(graph)的方法首先根据OCT图像中的区域项和边界项构造待分割的图，之后就可以将分割问题转换为图的最小割问题进行求解。Garvin等人提出一种结合图像边界和区域信息的三维图搜索(3D graph search)算法，其能够成功地分割出视网膜黄斑区域的5个边界。在该算法的基础上，Chiu等人结合动态规划理论对其进行了拓展，进一步提高了分割的准确率。虽然这些方法能够获得较高的准确率，但其算法复杂度较高，需要很长的运算时间。为了解决这些问题，Lee等人提出了一种可用于并行计算的图搜索分割算法，该方法在保证较高准确率的基础上，大大地缩短了算法运行的时间。Kafieh等人提出了一种基于扩散图(diffusion map)的分割算法，其能够分割出视网膜的12条边界。相比较其它基于图论的方法，该方法仅仅依赖于图像的区域信息，而不需要使用边界信息。而且，该算法在图像对比度和梯度较低的情况下野表现出较强的鲁棒性。上述基于三维图论的分割方法对噪声具有很强的鲁棒性，预处理步骤中对噪声的去除效果要求不高。但相比较其它类别的算法而言，大多数基于图论的算法的时间复杂度较高，同时在血管阴影较多的情况下，分割准确率有所下降。

发明内容

本发明提供了一种基于空间连续性约束的视盘区域OCT图像层次分割方法，其目的在于，为了克服视网膜视盘区域OCT图像中多个高反射区和血管阴影的存在对层次分割的影响，本发明能够准确分割视网膜中的ILM(Internal Limiting Membrane)、IS-OS(Inner segment-Outer segment)及BM(Bruch membrane)边界。

一种基于空间连续性约束的视盘区域OCT图像层次分割方法，包括以下步骤：

步骤1：对单只眼睛的各视网膜OCT图像进行降噪和ROI区域分割；

所述ROI区域是指视网膜OCT图像中的RNFL层和RPE层区域；

步骤2：将降噪后的视网膜OCT图像合成得到合成眼底图，从合成眼底图中定位血管阴影区域；

步骤3：使用A-scan分割算法依次对所述ROI区域进行分割，得到初步分割结果，利用血管阴影区域对每幅图像的初步分割结果进行修正；

将初步分割结果中位于血管阴影区域中的分割点剔除，再使用多项式拟合血管阴影区域两边的断点，得到连续的分割曲线，得到修正后的分割结果；

步骤4：按以下公式使用空间连续性约束优化步骤3得到的修正后的分割结果，获得ILM、IS-OS、BM的分割边界；

S＝{(i,j)||I_k(i,j)-I_k-1(i,j)|＜ε,gradient(i,j,k)＝1}

其中，I_k(i,j)表示第k幅OCT图像中位于第i行j列的初步分割结果点的坐标，slope(X,Y)为点X与点Y在图像坐标系上的斜率；集合S中为满足空间连续性约束条件的分割点；ε>0和ω>0分别表示两个相邻分割点之间欧式距离和梯度的阈值，ε∈(0,10)；ω∈(0,5)。

进一步地，采用BM3D算法对各视网膜OCT图像进行降噪。

进一步地，采用基于主动轮廓的FCM方法对降噪后的视网膜OCT图像进行分类，从分类结果中分割出ROI区域，具体过程如下：

步骤1.1：求解视网膜OCT图像中非相似性指标的价值函数的最小值，所述价值函数如下：

其中，H和W分别表示视网膜OCT图像的高和宽；J(·)为所要求解的目标函数，m是隶属度的加权指数，m∈[1,∞)，U为隶属度矩阵，U(i,j,k)表示视网膜OCT图像中位于第i行j列的像素点属于类别k的概率，U(i,j,k)＝p((i,j)∈μ_k)∈(0,1)；c表示聚类中心总数，取值为3；

μ_k表示第k个类别的聚类中心；d((i,j),μ_k)表示位于第i行j列的像素点与第k个聚类中心像素点的灰度值之差的二范数，d((i,j),μ_k)＝||I(i,j),I(μ_k)||₂；

m用于调节相似性度量与隶属度之间的权重关系；

所得3个分类结果图分别为：仅含有背景的分类结果、不含有RNFL层和RPE层的分类结果以及仅含有RNFL层和RPE层的分类结果图；

步骤1.2：选择图像中所有像素和最小的一个图像作为M，M为仅含有RNFL层和RPE层的分类图：

其中，U(k)表示第k类对应的图像；

步骤1.3：以ILM边界向下选择固定的像素数的方法将M中的RNFL层和RPE层进行分离，获得只含有RNFL层的图像M₁和只含有RPE层的图像M₂；

步骤1.4：采用主动轮廓模型对图像M₁和图像M₂进行细化处理，得到以二值图像C₁和二值图像C₂，二值图像C₁和二值图像C₂中的白色区域分别对应原始图像中RNFL层和RPE层的位置，基于二值图像C₁和二值图像C₂从原始图像中提取所需的RNFL层和RPE层构成的ROI区域。

进一步地，所述将降噪后的视网膜OCT图像合成得到合成眼底图的过程是指将所有降噪后的视网膜OCT图像在z轴方向上进行投影，以z轴方向灰度值的平均值作为合成图像中的灰度值。

将所有降噪后的视网膜OCT图像按列求解灰度值均值，该均值作为合成图像中的灰度值。例如，一张885*512的图像，经过按列求均值后变为1*512，共128张，则合成图像为128*512；

进一步地，所述从合成眼底图中定位血管阴影区域的过程如下：

对合成眼底图依次进行插值和滤波处理，再根据合成眼底图中血管区域特征，使用高斯滤波器对滤波处理后的合成眼底图像进行处理获得响应图像，最后使用大大津法对响应图像进行阈值分割，得到血管网络图；以OCT扫描线与血管网络图中的交点确定血管阴影区域。

响应图像中的像素亮度表示该像素属于血管的概率；

进一步地，使用A-scan分割算法依次对所述ROI区域进行分割的过程是指按照以下公式使用A-scan分割算法对所述二值化图像C₁、二值化图像C₂进行分割，得到二次分割结果：

其中，candidate(j)表示二值化图像中第j个列A-scan中属于边界的像素点的行序号，该像素点使目标函数达到最大值；g(i,j)表示二值化图像中位于第i行j列的像素点的梯度，distance(x,y)＝1/|x-y|表示两像素点x和y的行序号差值的倒数；α和γ分别为梯度项和距离项的系数，α∈[0,1]；γ∈[0,1]；

表示以i为变量，使f(i)达到最大值时，返回对应的i值的函数，f(i)＝α·g(i,j)+γ·distance(i,candidate(j-1))。

candidate(j)是初步分割结果，其中j表示图像中的列序号，candidate(j)表示图像中第j列的行序号。行序号和列序号确定一个点，该点即为初步分割的结果点，多个点组合成一条线。例如885*512的图像对应的初步分割结果为512个点，candidate(10)＝245,表示原始图像中点(245,10)为一个分割结果点。

有益效果

本发明提供了一种基于空间连续性约束的视盘区域OCT图像层次分割方法，该方法首先使用BM3D(Block-matching and 3D filtering)算法去除OCT图像中的散斑噪声；之后使用基于模糊C均值(FCM)和主动轮廓(C-V)的ROI(Region Of Interest)分割方法将图像中多个高反射区相互分离；然后使用ROI区域图像对血管阴影进行定位；针对原有A-scan分割算法易受到图像中血管阴影的影响，同时未考虑多张图像之间边界的连续性导致分割边界发生抖动的情况；然后使用ROI区域图像对血管阴影进行定位；使用A-scan分割算法依次对所述ROI区域进行分割，利用血管阴影区域对每幅图像的初步分割结果进行修正；最后使用空间连续性约束优化修正后的分割结果，获得ILM、IS-OS、BM的分割边界；应用本发明所述方法得到的最终分割结果与人工标定的金标准相比的有符号与无符号误差分别为：2.62±5.66(Mean±SD)μm和5.71±4.33(Mean±SD)μm。实验结果表明，该方法是一种行之有效的视盘区域OCT图像的层次分割方法，分割准确度较高，并对散斑噪声、血管阴影具有一定的鲁棒性。

附图说明

图1为本发明所述方法的流程图；

图2为FCM算法分割结果图，其中，(a)为原始图像，(b)为仅含有背景的分类结果图，(c)为不含有RNFL层和RPE层的分类结果图，(d)为仅含有RNFL层和RPE层的分类结果图；

图3为OCT图像血管阴影定位过程示意图，其中，(a)为OCT扫描图像，(b)为投影合成得到的眼底图，(c)为对合成眼底图进行血管分割的结果，上方横线为OCT系统扫描线，(d)为扫描线对应的OCT图；

图4为A-scan算法分割结果，其中，(a)为标注方框部分，(b)为(a)标注部分的图像放大细节图；

图5为血管阴影去除过程示意图，其中，(a)为血管阴影对分割结果造成的影响示意图，(b)为方框中标注的血管阴影的位置信息示意图，(c)为对血管阴影去除后的效果示意图；

图6为视网膜ONH区域OCT图像的三个边界分割示例，其中，(a)为第14个B-scan图像，(b)为第52个B-scan图像；

图7为视网膜ONH区域图像的分割示例，其中，(a)-(f)分别为第114、64、54、46、19、90个B-scan图像；

图8为应用本发明所述方法分割结果的三维重建展示，其中，(a)～(c)分别为ILM、ISOS以及BM边界，(d)为三个边界以3D的形式共同进行展示；

图9为本发明所述方法与专家标定的方法之间的Bland-Altman图展示，(a)～(c)分别为两种方法中ILM、ISOS、BM边界的分割结果对应的B&A图。

具体实施方式

下面将结合附图和实施例对本发明做进一步地说明。

本发明提出了一种基于空间连续性约束的视盘区域OCT图像层次分割方法，如图1所示，包括四个步骤：首先对图像进行降噪及ROI区域分割，以减弱散斑噪声和高反射区对分割结果的影响；之后对视网膜OCT图像中的血管阴影进行定位；然后结合血管阴影位置信息及A-scan分割算法对经过预处理后的图像进行分割，获得粗分割结果；最后使用空间连续性约束条件优化分割结果以获得最终更为准确的分割结果。

一种基于空间连续性约束的视盘区域OCT图像层次分割方法，具体包括以下步骤：

步骤1：对单只眼睛的各视网膜OCT图像进行降噪和ROI区域分割；

所述ROI区域是指视网膜OCT图像中的RNFL层和RPE层区域；

对OCT图像中的散斑噪声进行建模，如公式(1)：

v(x)＝u(x)+η(x),x∈Ω (1)

其中，u(x)表示理论上的无噪声图像，v(x)表示实际上我们看到的有噪声图像，η(x)表示噪声模型，用于描述图像中噪声的分布。可使用η(x)～N(0,σ²)的高斯噪声模型来近似代替图像中的噪声，然后使用BM3D算法进行降噪。

OCT图像中两个高反射区分别为RNFL层和RPE层。

采用基于主动轮廓的FCM方法对降噪后的视网膜OCT图像进行分类，从分类结果中分割出ROI区域，具体过程如下：

步骤1.1：分割为c类的过程等效为图像中非相似性指标的价值函数求解极小值的过程，其价值函数如公式(2)，求解视网膜OCT图像中非相似性指标的价值函数的最小值，所述价值函数如下：

其中，H和W分别表示视网膜OCT图像的高和宽；J(·)为所要求解的目标函数，m是隶属度的加权指数，m∈[1,∞)，U为隶属度矩阵，U(i,j,k)表示视网膜OCT图像中位于第i行j列的像素点属于类别k的概率，U(i,j,k)＝p((i,j)∈μ_k)∈(0,1)；c表示聚类中心总数，取值为3；

μ_k表示第k个类别的聚类中心；d((i,j),μ_k)表示位于第i行j列的像素点与第k个聚类中心像素点的灰度值之差的二范数，d((i,j),μ_k)＝||I(i,j),I(μ_k)||₂；

m用于调节相似性度量与隶属度之间的权重关系；

在本发明FCM分类过程中，设置c＝3。所得3个分类结果图分别为：仅含有背景的分类结果、不含有RNFL层和RPE层的分类结果以及仅含有RNFL层和RPE层的分类结果图，分别为图2中的(b)、(c)、(d)，选择图像中所有像素和最小的一个图像作为M，M即为仅含有RNFL层和RPE层的分类图(即图2(d))。

步骤1.2：选择图像中所有像素和最小的一个图像作为M，M为仅含有RNFL层和RPE层的分类图：

其中，U(k)表示第k类对应的图像；

步骤1.3：以ILM边界向下选择固定的像素数的方法将M中的RNFL层和RPE层进行分离，获得只含有RNFL层的图像M₁和只含有RPE层的图像M₂；

步骤1.4：采用主动轮廓模型对图像M₁和图像M₂进行细化处理，得到以二值图像C₁和二值图像C₂，二值图像C₁和二值图像C₂中的白色区域分别对应原始图像中RNFL层和RPE层的位置，基于二值图像C₁和二值图像C₂从原始图像中提取所需的RNFL层和RPE层构成的ROI区域。

步骤2：将降噪后的视网膜OCT图像合成得到合成眼底图，从合成眼底图中定位血管阴影区域；

所述将降噪后的视网膜OCT图像合成得到合成眼底图的过程是指将所有降噪后的视网膜OCT图像在z轴方向上进行投影，以z轴方向灰度值的平均值作为合成图像中的灰度值。

将所有降噪后的视网膜OCT图像按列求解灰度值均值，该均值作为合成图像中的灰度值。例如，一张885*512的图像，经过按列求均值后变为1*512，共128张，则合成图像为128*512；

如图3所示，所述从合成眼底图中定位血管阴影区域的过程如下：

对合成眼底图依次进行插值和滤波处理，再根据合成眼底图中血管区域特征，使用高斯滤波器对滤波处理后的合成眼底图像进行处理获得响应图像，如图3(b)，最后使用大大津法对响应图像进行阈值分割，得到血管网络图，如图3(c)；以OCT扫描线与血管网络图中的交点确定血管阴影区域。

响应图像中的像素亮度表示该像素属于血管的概率；

其中，使用的高斯滤波器定义如公式(4)：

其中，σ代表高斯分布的标准差。图3(c)中红线为OCT系统拍摄时扫描线的示意图，相干光扫过这条线得到图3(d)的图像。扫描线与血管网络的交点与OCT图像中的血管阴影一一对应，通过这种对应关系可以准确地确定图像中的阴影位置。

步骤3：使用A-scan分割算法依次对所述ROI区域进行分割，得到初步分割结果，利用血管阴影区域对每幅图像的初步分割结果进行修正；

将初步分割结果中位于血管阴影区域中的分割点剔除，再使用多项式拟合血管阴影区域两边的断点，得到连续的分割曲线，得到修正后的分割结果；

使用A-scan分割算法依次对所述ROI区域进行分割的过程是指按照以下公式使用A-scan分割算法对所述二值化图像C₁、二值化图像C₂进行分割，得到二次分割结果：

其中，candidate(j)表示二值化图像中第j个列A-scan中属于边界的像素点的行序号，该像素点使目标函数达到最大值；g(i,j)表示二值化图像中位于第i行j列的像素点的梯度，distance(x,y)＝1/|x-y|表示两像素点x和y的行序号差值的倒数；α和γ分别为梯度项和距离项的系数，α∈[0,1]；γ∈[0,1]；

表示以i为变量，使f(i)达到最大值时，返回对应的i值的函数，f(i)＝α·g(i,j)+γ·distance(i,candidate(j-1))。

通过对公式(5)求解，可以得到视网膜层次的初步分割结果。在该结果中，血管阴影的存在会导致分割边界误差较大且分割边界不平滑，如图4所示。

通过观察可以发现，血管阴影的对分割结果的影响多局限于阴影区域，在阴影区域以外则对结果无影响；因此，结合步骤2得到的血管阴影位置信息，将位于阴影区域的分割点剔除，然后平滑阴影区域两边的断点，得到连续的分割曲线，分割结果如图5所示。

candidate(j)是初步分割结果，其中j表示图像中的列序号，candidate(j)表示图像中第j列的行序号。行序号和列序号确定一个点，该点即为初步分割的结果点，多个点组合成一条线。例如885*512的图像对应的初步分割结果为512个点，candidate(10)＝245,表示原始图像中点(245,10)为一个分割结果点。

步骤4：按以下公式使用空间连续性约束优化步骤3得到的修正后的分割结果，获得ILM、IS-OS、BM的分割边界；

其中，I_k(i,j)表示第k幅OCT图像中位于第i行j列的初步分割结果点的坐标，slope(X,Y)为点X与点Y在图像坐标系上的斜率；集合S中为满足空间连续性约束条件的分割点；ε>0和ω>0分别表示两个相邻分割点之间欧式距离和梯度的阈值，ε∈(0,10)；ω∈(0,5)。优化后的分割结果如图6所示，在图6中，实线为本文方法的分割结果，虚线专家标定的金标准；边界1为ILM，边界2为IS-OS，边界3为BM。

本发明所提出的视盘区域OCT层次分割方法利用MATLAB R2013a实现以及对实验数据进行分析，计算误差、绘制Bland-Altman图。

为验证本发明提出的方法的性能，从实验数据集中随机选取20张图像进行分割并分析误差。其中实验数据全部使用Topcon 3D OCT-1Maestro系统进行采集，每只眼睛采集的图像数据维度为512*128*885，对应生理组织上的大小为6mm*6mm*2mm。

表1为三个分割边界与金标准之间的误差统计，其中ILM边界的分割结果最为准确，有符号和无符号误差分别为(-1.52±3.70um)和(3.36±2.47um)。第三个边界为BM边界，位于RPE层和脉络膜之间，由于脉络膜区域的组织反射率不均衡，导致该边界的分割误差较大。图7列出了应用本发明所述方法得到的分割结果与专家手动分割结果的对比图，每幅子图的上面一行是专家标定的结果，下面一行是本文方法的分割结果，图8对层次分割结果的三维重建示意图。从结果中可以看出，本发明所述的方法可以准确地分割出相应的视网膜边界。

表1本文方法分割结果与金标准之间的平均误差与标准差(Mean±SD，单位为um)

*统计误差时不包含Neural Canal区域的分割结果

表2为本发明方法与Lee等人和Saleh等人方法针对ILM边界分割结果的比较结果。由表中内容可知，相比较于其它算法，本发明分割结果有更小的有符号与无符号误差。然而同时也发现应用本发明所述方法得到的分割结果有更大的标准差，特别是在有符号误差统计中。经过实验分析发现，图像中ONH区域的ILM边界比较陡峭，在分割时相比较视网膜其它区域，会有较大的偏差。同时在统计误差的过程中，本发明所述方法采用了计算分割的边界与金标准之间垂直像素差的方式，这样的方法会引起统计结果方差的偏大。

表2边界分割的平均有符号与无符号误差(Mean±SD，单位为um)

*统计误差时包含Neural Canal区域的分割部分

为了评估本发明方法的分割结果与金标准之间的一致性，使用Bland-Altman(B&A)图对两种方法的结果进行分析，如图9所示。其中，Y轴为两种方法对同一位置的测量值的差值，X轴为两种测量值的均值，图中Mean和SD分别为整体测量值差值的均值和标准差。B&A图中，如果有95％以上的数据点落在±1.96SD区域内，则可以认为两种方法具有一致性，即两种方法在测量时可以交换使用。由图9可知，三个分割边界对应的B&A图中分别有96.2％、94.7％和99.7％的点落在该区域。因此，通过对结果的分析可知，本发明方法在ONH区域视网膜层次分割上可以替代专家标定的方法，在很大程度上提高分层的效率以及减少手工分割的工作量。

本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。

Claims (6)

1.一种基于空间连续性约束的视盘区域OCT图像层次分割方法，其特征在于，包括以下步骤：

步骤1：对单只眼睛的各视网膜OCT图像进行降噪和ROI区域分割；

所述ROI区域是指视网膜OCT图像中的RNFL层和RPE层区域；

步骤2：将降噪后的视网膜OCT图像合成得到合成眼底图，从合成眼底图中定位血管阴影区域；

步骤3：使用A-scan分割算法依次对所述ROI区域进行分割，得到初步分割结果，利用血管阴影区域对每幅图像的初步分割结果进行修正；

将初步分割结果中位于血管阴影区域中的分割点剔除，再使用多项式拟合血管阴影区域两边的断点，得到连续的分割曲线，得到修正后的分割结果；

步骤4：按以下公式使用空间连续性约束优化步骤3得到的修正后的分割结果，获得ILM、IS-OS、BM的分割边界；

S＝{(i,j)||I_k(i,j)-I_k-1(i,j)|＜ε,gradient(i,j,k)＝1}

其中，I_k(i,j)表示第k幅OCT图像中位于第i行j列的初步分割结果点的坐标，slope(X,Y)为点X与点Y在图像坐标系上的斜率；集合S中为满足空间连续性约束条件的分割点；ε>0和ω>0分别表示两个相邻分割点之间欧式距离和梯度的阈值，ε∈(0,10)；ω∈(0,5)。

2.根据权利要求1所述的方法，其特征在于，采用BM3D算法对各视网膜OCT图像进行降噪。

3.根据权利要求2所述的方法，其特征在于，采用基于主动轮廓的FCM方法对降噪后的视网膜OCT图像进行分类，从分类结果中分割出ROI区域，具体过程如下：

步骤1.1：求解视网膜OCT图像中非相似性指标的价值函数的最小值，所述价值函数如下：

其中，H和W分别表示视网膜OCT图像的高和宽；J(·)为所要求解的目标函数，m是隶属度的加权指数，m∈[1,∞)，U为隶属度矩阵，U(i,j,k)表示视网膜OCT图像中位于第i行j列的像素点属于类别k的概率，U(i,j,k)＝p((i,j)∈μ_k)∈(0,1)；c表示聚类中心总数，取值为3；

μ_k表示第k个类别的聚类中心；d((i,j),μ_k)表示位于第i行j列的像素点与第k个聚类中心像素点的灰度值之差的二范数，d((i,j),μ_k)＝||I(i,j),I(μ_k)||₂；

步骤1.2：选择图像中所有像素和最小的一个图像作为M，M为仅含有RNFL层和RPE层的分类图：

其中，U(k)表示第k类对应的图像；

步骤1.3：以ILM边界向下选择固定的像素数的方法将M中的RNFL层和RPE层进行分离，获得只含有RNFL层的图像M₁和只含有RPE层的图像M₂；

步骤1.4：采用主动轮廓模型对图像M₁和图像M₂进行细化处理，得到以二值图像C₁和二值图像C₂，二值图像C₁和二值图像C₂中的白色区域分别对应原始图像中RNFL层和RPE层的位置，基于二值图像C₁和二值图像C₂从原始图像中提取所需的RNFL层和RPE层构成的ROI区域。

4.根据权利要求2所述的方法，其特征在于，所述将降噪后的视网膜OCT图像合成得到合成眼底图的过程是指将所有降噪后的视网膜OCT图像在z轴方向上进行投影，以z轴方向灰度值的平均值作为合成图像中的灰度值。

5.根据权利要求4所述的方法，其特征在于，所述从合成眼底图中定位血管阴影区域的过程如下：

对合成眼底图依次进行插值和滤波处理，再根据合成眼底图中血管区域特征，使用高斯滤波器对滤波处理后的合成眼底图像进行处理获得响应图像，最后使用大大津法对响应图像进行阈值分割，得到血管网络图；以OCT扫描线与血管网络图中的交点确定血管阴影区域。

6.根据权利要求3所述的方法，其特征在于，使用A-scan分割算法依次对所述ROI区域进行分割的过程是指按照以下公式使用A-scan分割算法对所述二值化图像C₁、二值化图像C₂进行分割，得到二次分割结果：

其中，candidate(j)表示二值化图像中第j个列A-scan中属于边界的像素点的行序号，该像素点使目标函数达到最大值；g(i,j)表示二值化图像中位于第i行j列的像素点的梯度，distance(x,y)＝1/|x-y|表示两像素点x和y的行序号差值的倒数；α和γ分别为梯度项和距离项的系数，α∈[0,1]；γ∈[0,1]；

表示以i为变量，使f(i)达到最大值时，返回对应的i值的函数，f(i)＝α·g(i,j)+γ·distance(i,candidate(j-1))。

Images