CN108269258B - 从oct角膜图像分割角膜结构的方法和*** - Google Patents

Abstract

本发明提供了一种从OCT角膜图像分割角膜结构的方法和***，涉及生物医学图像技术领域，该方法包括：对OCT角膜图像进行缩放处理；确定角膜顶点的横坐标和角膜各层顶点的初始纵坐标；以圆拟合方式确定初始分界面；按照设定的邻域宽度对角膜各层的初始分界面进行边界扩充，得到环形感兴趣区；得到矩形感兴趣区；对矩形感兴趣区分别进行平均灰度差分处理和扁平化处理的逆处理，得到修正分界面；基于垂直平分线过圆心原理，对每个修正分界面分别进行霍夫变换圆拟合操作，得到角膜各层的圆心和半径；分别对角膜各层的修正分界面进行平抛运动的卡尔曼滤波。本发明提高了OCT角膜图像分割角膜各层的速度和精度。

Description

从OCT角膜图像分割角膜结构的方法和***

技术领域

本发明涉及生物医学图像技术领域，尤其是涉及一种从OCT角膜图像分割角膜结构的方法和***。

背景技术

光学相干断层成像(Optical Coherence Tomography,OCT)是一种对成像表面组织结构无损的光学信号获取和处理方法。近些年来，OCT已经成为一种重要的具有创新意义的角膜测量和诊断技术。由于对角膜的不同层进行量化分析可以帮助诊断视功能改变并帮助医生对角膜水肿、眼球高压和内皮功能缺陷进行诊断，以及屈光手术术前、术后的角膜厚度分析。因此，在保证速度和精度的前提下，对角膜的不同层进行分割，获得厚度和曲率信息具有重要的意义。

然而，现有的全自动角膜分割方法存下以下问题：

对眼前段OCT产生的低信噪比区域，如角膜中央过饱和伪影、水平方向伪影或者由于远心扫描模式(目前大多数频域OCT采用该方式)引起的角膜两侧信号弱，没有进行有效处理；现有方法速度普遍较慢，影响医生后续的实时诊断，无法同时满足临床的实时性和精确性需求。

发明内容

有鉴于此，本发明的目的在于提供一种从OCT角膜图像分割角膜结构的方法和***，以提高OCT角膜图像分割角膜各层的速度和精度。

第一方面，本发明实施例提供了一种从OCT角膜图像分割角膜结构的方法，包括：根据OCT角膜图像的OCT设备的分辨率，对OCT角膜图像进行缩放处理，得到横向分辨率和纵向分辨率相同的OCT等分辨率图像；其中，所述OCT角膜图像为所述OCT设备生成的二维灰度图；根据所述OCT等分辨率图像中各像素的灰度值确定角膜顶点的横坐标和角膜各层顶点的初始纵坐标，其中，所述角膜各层顶点包括：上皮层顶点、鲍尔曼层顶点和内皮层顶点；根据所述角膜顶点的横坐标和所述角膜各层顶点的初始纵坐标，以圆拟合方式确定所述角膜各层的初始分界面；按照设定的邻域宽度对所述角膜各层的初始分界面进行边界扩充，得到所述角膜各层的环形感兴趣区；对所述角膜各层的环形感兴趣区分别进行扁平化处理，得到所述角膜各层的矩形感兴趣区；对所述角膜各层的矩形感兴趣区分别进行平均灰度差分处理和所述扁平化处理的逆处理，得到所述角膜各层的修正分界面；对于所述角膜各层的修正界面中的每个修正分界面，基于垂直平分线过圆心原理，对每个所述修正分界面分别进行霍夫变换圆拟合操作，得到所述角膜各层的圆心和半径；基于所述角膜各层的圆心和半径，分别对所述角膜各层的修正分界面进行平抛运动的卡尔曼滤波，得到所述角膜各层的优化分界面。

进一步地，根据OCT角膜图像的OCT设备的分辨率，对OCT角膜图像进行缩放处理，得到横向分辨率和纵向分辨率相同的OCT等分辨率图像包括：获取OCT角膜图像的OCT设备的分辨率，其中，所述OCT设备的横向分辨率记为A1，轴向分辨率记为B1；获取所述OCT角膜图像的列数和行数，所述列数记为C1，所述行数记为D1；设置OCT等分辨率图像的列数为：C1×A1/B1，行数为D1；或者设置OCT等分辨率图像的列数为：C1，行数为D1×B1/A1；以使所述OCT等分辨率图像的横向分辨率和纵向分辨率相同。

进一步地，根据所述OCT等分辨率图像中各像素的灰度值确定角膜顶点的横坐标和角膜各层顶点的初始纵坐标包括：对所述OCT等分辨率图像分别进行轴向投影和横向投影；根据所述轴向投影得到的灰度值大小确定角膜顶点的横坐标；根据所述轴向投影得到的灰度值大小确定所述OCT等分辨率图像的中央伪影区；将所述横向投影得到的最大灰度值对应区域，确定为所述OCT等分辨率图像的水平伪影区；基于确定出的中央伪影区和水平伪影区，获取角膜各层顶点的初始纵坐标。

进一步地，根据所述轴向投影得到的灰度值大小确定所述OCT等分辨率图像的中央伪影区包括：将所述OCT等分辨率图像沿横向坐标轴等宽度地分为三个区域，依次为区域I、区域Ⅱ和区域Ⅲ，根据所述轴向投影计算区域I和Ⅲ的平均灰度，分别记为m_I和m_III，并设定阈值T_cent为：

其中，g_max为所述轴向投影中的最大值，所述最大值对应的横坐标为P_max，ρ为预先设定的值；对所述轴向投影得到的灰度值，以P_max对应的位置为中心，分别向两侧逐一查找小于T_cent的灰度值，将两侧中首个查找到的灰度值对应的位置分别记为p_left和p_right；将所述OCT等分辨率图像中p_left和p_right之间的对应区域确定为中央伪影区。

进一步地，基于确定出的中央伪影区和水平伪影区，获取角膜各层顶点的初始纵坐标包括：对所述横向投影组成的向量进行一阶前向差分运算，将运算结果中的最大向量对应的纵坐标记为所述角膜的上皮层顶点的纵坐标；将所述OCT等分辨率图像中的中央伪影区的每列像素的灰度值设置为0，对所述上皮层顶点以下区域进行扁平化操作，得到所述OCT等分辨率图像的局部拉伸图像；对所述局部拉伸图像进行再次横向投影，对所述再次横向投影组成的向量进行一阶前向差分运算，将运算结果中的最大向量对应的纵坐标记为所述角膜的鲍尔曼层顶点的纵坐标；将运算结果中的最小向量对应的纵坐标记为所述角膜的内皮层顶点的纵坐标。

进一步地，确定出所述OCT等分辨率图像的水平伪影区之后，所述方法还包括：按照以下公式对所述水平伪影区进行抑制操作：

其中，f(x,y)为所述OCT等分辨率图像中像素(x,y)处的灰度值，W为所述OCT等分辨率图像的横向长度、H为所述OCT等分辨率图像的纵向长度，x与y为非负整数，取值范围分别为[0,W-1]与[0,H-1]其中min与max分别为

的最小值与最大值。

进一步地，对所述角膜各层的矩形感兴趣区分别进行平均灰度差分处理包括：采用以下公式分别计算所述角膜各层的矩形感兴趣区的界面：

其中，j的取值范围为[1,H_ROI-1]，H_ROI为扁平化后矩形感兴趣区的高度；f(x,y)为所述矩形感兴趣区中像素(x,y)处的灰度值；所述角膜各层的矩形感兴趣区包括：空气-上皮分界面对应的矩形感兴趣区、上皮层-鲍尔曼层分界面对应的矩形感兴趣区和内皮层-房水分界面对应的矩形感兴趣区；将空气-上皮分界面对应的矩形感兴趣区计算出的每列D_j中的最大值作为上皮层界面在该列的纵坐标；将上皮层-鲍尔曼层分界面对应的矩形感兴趣区计算出的每列D_j中的最大值作为鲍尔曼层界面在该列的纵坐标；将内皮层-房水分界面对应的矩形感兴趣区计算出的每列D_j中的最小值作为内皮层界面在该列的纵坐标。

进一步地，基于垂直平分线过圆心原理，对每个所述修正分界面分别进行霍夫变换圆拟合操作，得到所述角膜各层的圆心和半径包括：对所述角膜各层的每个修正分界面，均进行如下操作得到所述角膜各层的圆心和半径：按照设定的第一间隔，依次从当前修正分界面上选取三个像素点K、M和N，连接KM和MN；将KM和MN各自垂直平分线的交点确定为当前圆心；按顺序选取K、M和N，以及圆心的确定过程，直至所述当前修正分界面上各像素点均被选取过为止，将重合次数最多的圆心作为所述当前修正分界面的最终圆心；以所述最终圆心为起点，按照设定的第二间隔，从所述当前修正分界面上逐一选取像素点，计算选取的像素点到所述最终圆心的距离，将所述距离作为当前半径；将最多重复次数的半径确定为所述当前修正分界面的半径。

进一步地，基于所述角膜各层的圆心和半径，分别对所述角膜各层的修正分界面进行平抛运动的卡尔曼滤波，得到所述角膜各层的优化分界面，包括：分别以所述角膜各层的顶点为原点建立坐标系，设置下述抛物线方程以获取所述角膜各层的优化分界面：

y＝ax²/2

其中，

a为抛物线方程的系数，该系数还用于后续确定卡尔曼滤波模型的状态方程，x_N为各个所述修正分界面对应的圆与所述OCT等分辨率图像左边界或右边界交点的横向坐标，r为所述修正分界面对应的圆的半径；按照如下方式设置卡尔曼滤波模型中的状态参数Ф和方差P：

φ(k)＝Aφ(k-1)+Ba

其中，

是二维向量，

y是当前角膜层的优化分界面上各点的纵坐标值，v_y是相邻边界点的速度，下标k表示当前边界点，下标(k-1)表示上一个边界点；A是状态转移矩阵，B是控制输入矩阵，^表示预测值，T表示转置；

式中dt表示从前一边界点至当前边界点的空间(时间)间隔；构建卡尔曼滤波模型中的校正方程如下：

其中，H＝[1 0]，R是预先设定的测量噪声的标准差，K将随着R和

进行更新；设置状态参数Ф按如下公式进行更新：

其中，Z(k)为所述当前角膜层的修正分界面上的k点的纵坐标值；设置方差P按照如下公式更新：

将更新后的Ф和P代入上述校正方程，直至方差P趋近于最小方差值。

第二方面，本发明实施例还提供一种从OCT角膜图像分割角膜结构的***，包括：分辨率调整模块，用于根据OCT角膜图像的OCT设备的分辨率，对OCT角膜图像进行缩放处理，得到横向分辨率和纵向分辨率相同的OCT等分辨率图像；其中，所述OCT角膜图像为所述OCT设备生成的二维灰度图；顶点坐标确定模块，用于根据所述OCT等分辨率图像中各像素的灰度值确定角膜顶点的横坐标和角膜各层顶点的初始纵坐标，其中，所述角膜各层顶点包括：上皮层顶点、鲍尔曼层顶点和内皮层顶点；初始界面确定模块，用于根据所述角膜顶点的横坐标和所述角膜各层顶点的初始纵坐标，以圆拟合方式确定所述角膜各层的初始分界面；环形感兴趣区确定模块，用于按照设定的邻域宽度对所述角膜各层的初始分界面进行边界扩充，得到所述角膜各层的环形感兴趣区；矩形感兴趣区确定模块，用于对所述角膜各层的环形感兴趣区分别进行扁平化处理，得到所述角膜各层的矩形感兴趣区；修正分界面确定模块，用于对所述角膜各层的矩形感兴趣区分别进行平均灰度差分处理和所述扁平化处理的逆处理，得到所述角膜各层的修正分界面；霍夫变换圆拟合模块，用于对于所述角膜各层的修正界面中的每个修正分界面，基于垂直平分线过圆心原理，对每个所述修正分界面分别进行霍夫变换圆拟合操作，得到所述角膜各层的圆心和半径；卡尔曼滤波模块，用于基于所述角膜各层的圆心和半径，分别对所述角膜各层的修正分界面进行平抛运动的卡尔曼滤波，得到所述角膜各层的优化分界面。

本发明实施例提供的从OCT角膜图像分割角膜结构的方法和***，首先对角膜图像进行缩放及坐标确定，再通过平均灰度差分方法确定角膜各层分界面，再通过抛物线内切圆确定抛物线系数，进而采用平抛运动建模的卡尔曼滤波来对平均灰度差分方法获得的分界面进行修正和约束，可以快速并精确地获得角膜图像各层的边界，提高了OCT角膜图像分割角膜各层的速度和精度。

为使本发明的上述目的、特征和优点能更明显易懂，下文特举较佳实施例，并配合所附附图，作详细说明如下。

附图说明

为了更清楚地说明本发明具体实施方式或现有技术中的技术方案，下面将对具体实施方式或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图是本发明的一些实施方式，对于本领域技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明实施例一提供的从OCT角膜图像分割角膜结构的方法的流程示意图；

图2(a)为本发明实施例一提供的角膜图像；

图2(b)为图2(a)提供的角膜图像的轴向投影；

图3(a)为本发明实施例一提供的角膜图像；

图3(b)为图3(a)提供的角膜图像的横向投影；

图4(a)为本发明实施例一提供的角膜图像；

图4(b)为图4(a)提供的角膜图像减少中央伪影区影响的示意图；

图4(c)为图4(b)提供的示意图扁平化后的示意图；

图4(d)为图4(a)提供的角膜图像的横向投影；

图5(a)为本发明实施例一提供的空气-上皮分界面ROI示意图；

图5(b)为本发明实施例一提供的空气-上皮分界面ROI扁平化示意图；

图6(a)为本发明实施例一提供的采用平均灰度差分计算空气-上皮分界面并对其进行Savitzky-Golay滤波后的示意图；

图6(b)为图6(a)提供的示意图进行逆扁平化后的示意图；

图7(a)为本发明实施例一提供的垂直平分弦的线段相交于圆心的示意图

图7(b)为本发明实施例一提供的用于对圆心进行投票的序列等间距三元组的示意图，

图8(a)为本发明实施例一提供的霍夫变换的圆心投票直方图；

图8(b)为本发明实施例一提供的霍夫变换的半径投票直方图；

图9为本发明实施例一提供的霍夫变换获得的空气-上皮分界面示意图；

图10(a)图10(b)为本发明实施例一提供的平抛运动示意图；

图11为本发明实施例一提供的角膜各层的优化分界面示意图；

图12为本发明实施例二提供的从OCT角膜图像分割角膜结构的***的示意图。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合附图对本发明的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

下面结合附图，对本发明的具体实施方式作详细说明。

实施例一

本发明实施例一提供了一种从OCT角膜图像分割角膜结构的方法，参见图1所示的从OCT角膜图像分割角膜结构的方法的流程示意图，该方法包括如下步骤：

步骤S11，根据OCT角膜图像的OCT设备的分辨率，对OCT角膜图像进行缩放处理，得到横向分辨率和纵向分辨率相同的OCT等分辨率图像。

其中，OCT角膜图像为OCT设备生成的二维灰度图。由于不同的OCT成像设备在横向和轴向的分辨率不同，考虑到后续可对中央角膜进行圆拟合，图像需要进行缩放以使图像分辨率在横轴与纵轴相同，从而使得圆在新的图像上依旧为圆。缩放过程可以保持图像高度不变，通过缩放宽度使缩放后的图像横向和纵向像素分辨率相同；也可以保持图像宽度不变，通过缩放高度进行；具体的缩放步骤如下：

获取OCT角膜图像的OCT设备的分辨率，其中OCT设备的横向分辨率记为A1，轴向分辨率记为B1；获取所述OCT角膜图像的列数和行数，列数记为C1，行数记为D1；设置OCT等分辨率图像的列数为：C1×A1/B1，行数为D1；或者设置OCT等分辨率图像的列数为：C1，行数为D1×B1/A1；以使所述OCT等分辨率图像的横向分辨率和纵向分辨率相同。

在后续步骤对缩放图像进行分割之后，图像可逆缩放至原始尺寸，以得到原图像的分割结果。

步骤S12，根据OCT等分辨率图像中各像素的灰度值确定角膜顶点的横坐标和角膜各层顶点的初始纵坐标。

具体地，可以对等分辨率图像进行轴向投影和横向投影，根据投影图像的灰度分布确定上述横坐标和纵坐标。其中，上述角膜各层顶点包括：上皮层顶点、鲍尔曼层顶点和内皮层顶点。正常人眼球是接近球状的。按照Gullstrand的模拟眼模型，人眼沿光轴的角膜切面是两个圆，OCT采用的是沿着眼主光轴的切面成像方法，因此上述角膜各层顶点的横坐标均使用角膜顶点的横坐标。

步骤S13，根据上述角膜顶点的横坐标和角膜各层顶点的初始纵坐标，以圆拟合方式确定角膜各层的初始分界面。

由于正常人眼角膜为圆形，根据已得到的横纵坐标和已知的角膜切面的半径可使用圆拟合方式确定角膜各层的初始分界面。

步骤S14，按照设定的邻域宽度对角膜各层的初始分界面进行边界扩充，得到角膜各层的环形感兴趣区。

上述邻域宽度指沿纵轴方向的宽度，为了保证将层边界包括至环形感兴趣区的同时环形感兴趣区不超出暗到亮或亮到暗区域，并避免增加计算复杂度，针对角膜各层可以分别设定合适的邻域宽度，该邻域宽度可以是定值或者在一定区间内变化的值，例如：15像素-35像素之间的数值，优选30像素的邻域宽度。

步骤S15，对角膜各层的环形感兴趣区分别进行扁平化处理，得到角膜各层的矩形感兴趣区。

采用扁平化操作的目的是让后续的平均灰度差分方法获得的边界更准确。扁平化操作以各分界面的环形感兴趣区的中心线和顶点纵坐标为基准，将其中的每一列像素向上移动，得到水平的矩形感兴趣区。

步骤S16，对角膜各层的矩形感兴趣区分别进行平均灰度差分处理和上述扁平化处理的逆处理，得到角膜各层的修正分界面。

步骤S17，对于角膜各层的修正界面中的每个修正分界面，基于垂直平分线过圆心原理，对每个修正分界面分别进行霍夫变换圆拟合操作，得到角膜各层的圆心和半径。

霍夫变换是一种重要的圆检测方法，在本实施例中，采用弦垂直平分线通过圆心原理的霍夫变换来加速投票过程，使得算法的复杂度降低。

步骤S18，基于上述角膜各层的圆心和半径，分别对角膜各层的修正分界面进行平抛运动的卡尔曼滤波，得到角膜各层的优化分界面。

具体地，通过基于平抛运动的卡尔曼滤波对步骤S16中得到的角膜各层的修正分界面进行修正和约束，其中使用步骤S17中霍夫变换获得的圆作为抛物线的内切圆，由该圆的圆心和半径确定上述卡尔曼滤波使用的抛物线的参数。

本发明实施例提供的从OCT角膜图像分割角膜结构的方法，首先对角膜图像进行缩放及坐标确定，再通过平均灰度差分确定角膜各层分界面，再通过抛物线内切圆确定抛物线系数，进而采用平抛运动建模的卡尔曼滤波来对平均灰度差分方法获得的分界面进行修正和约束，可以快速并精确地获得角膜图像各层的边界，，提高了OCT角膜图像分割角膜各层的速度和精度。

为了在保持精度的同时降低计算复杂度，本实施例先对图像做预处理，其中角膜图像如图2(a)所示。预处理分为四个步骤：图像缩放，定位角膜顶点的横坐标，定位角膜各层顶点的纵坐标和水平噪声的抑制。

一、图像缩放

由于不同的OCT成像设备在横向和轴向的分辨率不同，考虑到后续可对中央角膜进行圆拟合，图像需要进行缩放以使像素大小在横轴与纵轴相同，从而使得圆在新的图像上依旧为圆。缩放过程可以保持图像高度不变，通过缩放宽度使缩放后的图像横向和纵向像素分辨率相同；也可以保持图像宽度不变，通过缩放高度进行。

例如，在角膜图像集中，每幅图像包括1000列，每列的像素个数为512个，横向和轴向的组织采样分辨率分别为6.1μm和4.6μm。为了便于圆检测，我们首先将图像放大到1326(1000x6.1/4.6)(横向)像素x 512(轴向)像素使得横向和纵向的像素分辨率相同。在后续步骤对缩放图像进行分割之后，图像可逆缩放至原始尺寸，以得到原图像的分割结果。

二、通过轴向投影来定位角膜顶点的横坐标

具体包括：

1.对图像沿轴向进行投影，投影公式如下：

轴向投影后的结果如图2(b)所示。

2.计算轴向投影的最大值g_max以及对应的横坐标p_max。

3.将图2(a)中的图像沿横向坐标轴等宽度地分为三个区域，依次为区域I、区域Ⅱ和区域Ⅲ，并假定中央伪影不会出现在区域I和Ⅲ中。

4.根据上述轴向投影计算区域I和Ⅲ的平均灰度，分别记为m_I和m_III，并设定阈值T_cent为：

其中，g_max为轴向投影的最大值，该最大值对应的横坐标为P_max，ρ为预先设定的值，在本实施例中ρ的取值为0.5左右。

5.以P_max对应的位置为中心，分别向两侧逐一查找小于T_cent的灰度值，将两侧中首个查找到的灰度值对应的位置分别记为p_left和p_right，将p_left和p_right之间的对应区域确定为中央伪影区。

6.计算角膜顶点的横坐标A_x：

A_x＝(p_left+p_right)/2 (3)

参见图2(a)和图2(b)所示，其中三条直线从左到右依次为中央伪影区的左边界、角膜顶点所在直线和中央伪影区的右边界。

三、通过横向投影来定位角膜各层顶点的纵坐标

其中，定位上皮层顶点纵坐标具体包括：

1.参见图3(a)所示的角膜图像，对其进行横向投影，投影公式如下：

得到的投影曲线如图3(b)中上侧的曲线所示，需要注意的是图3(b)的横坐标轴对应于图3(a)的纵坐标轴。由图可知，投影曲线在空气-上皮分界面顶点处获得极大值(上皮层顶点处)，将横向投影得到的最大值对应区域(水平伪影处)，确定为水平伪影区。

2.由于水平伪影的影响，投影曲线在空气-上皮分界面顶点处不一定是最大值，为了更加精确地获得角膜顶点，可对上述投影公式获得的投影向量作一阶前向差分运算，公式如下：

g'(y)＝g(y+1)-g(y) (5)

横向投影向量一阶前向差分运算结果如图3(b)中下侧的曲线所示，其中差分曲线在空气-上皮分界面顶点处获得最大值，对应的y坐标记为

对应于上皮层顶点的纵坐标。

定位鲍尔曼层顶点和内皮层顶点的纵坐标，具体包括：

1.参见图4(a)所示的角膜图像，为了减少中央伪影的影响，将中央伪影区的每列像素的灰度值设置为0，如图4(b)所示。

2.对上皮层顶点以下区域进行扁平化操作，得到局部拉伸图像，如图4(c)所示，具体采取如下方式：对每一列中像素，将其向上移动C_epithe(y)-A′_y个位置，其中A′_y为空气-上皮层分界面顶点的纵坐标，C_epithe(y)为空气-上皮层界面的纵坐标；扁平化操作之后，空气-上皮层变为水平。

3.对扁平化后的图像按照公式(4)进行横向投影，如图4(d)中上侧的曲线所示，需要注意的是图4(d)的横坐标轴对应于图4(a)的纵坐标轴。

4.按照公式(5)作前向差分运算，如图4(d)中下侧的曲线所示，将运算结果中的最大向量对应的纵坐标记为鲍尔曼层顶点的纵坐标；将运算结果中的最小向量对应的纵坐标记为内皮层顶点的纵坐标。

四、水平噪声的抑制

在对图像预处理过程中，水平噪声的抑制指对水平伪影的抑制操作，具体可通过如下公式进行：

其中，f(x,y)为图像中像素(x,y)处的灰度值，W为图像的横向长度、H为图像的纵向长度，x与y为非负整数，取值范围分别为[0,W-1]与[0,H-1]其中min与max分别为

的最小值与最大值。

上述公式本质上是将每点的灰度变为该点的灰度减去该行灰度的平均然后线性拉伸到0-255灰度，这样确保那些灰度类似的行(对应于水平伪影)变成接近0，从而消除水平伪影。

在完成对角膜图像的预处理后，采用平均灰度差分提取各层修正分界面。由于角膜各层分界面具有层状边界结构，且在深度方向存在着亮到暗或者暗到亮的灰度过渡，因此可利用这两个特性来帮助分割角膜各层。具体可通过三个步骤实现：分界面的感兴趣区(Region of Interest，ROI)提取，感兴趣区域扁平化和采用平均灰度差分方法获取修正分界面，具体如下：

一、分界面的感兴趣区域提取

为了使后续的边界提取更准确，本实施例采用了对每层边界邻域即ROI进行提取和扁平化。ROI提取需要确定ROI的中心线位置、ROI的宽度和高度方向邻域半径。由各分界面顶点坐标(A_x,A_y)和半径可以确定唯一的圆C_init，该圆每列的纵坐标为C_init(y)。以空气-上皮分界面为例，由于前文预处理过程已经估算出空气-上皮分界面顶点的横坐标(即角膜顶点的横坐标)A_x和纵坐标

根据先验知识又已知空气-上皮分界面所形成的圆半径约为7.7mm，因此可以通过圆拟合方式得到一个圆

该圆每列的纵坐标为

该圆大致接近角膜上皮层，

即为空气-上皮分界面ROI的中心线。对其他两层分界面ROI的中心线的提取可以采用类似的方法，横坐标为A_x，纵坐标提取见前文预处理过程，根据先验知识,内皮-房水分界面所形成的圆半径约为6.8mm，上皮-鲍尔曼分界面所形成的圆半径约为7.6mm。ROI的宽度与缩放后的图像宽度W相同，ROI的高度方向邻域半径需根据不同边缘类型设定。邻域半径太大，会增加计算复杂度，并且可能超出了暗到亮或亮到暗区域的范围，导致后续的平均灰度差分提取分界面精度下降。由于根据各分界面顶点坐标和先验半径确定的初始圆有可能存在误差，邻域半径太小可能使得层分界面没有包括在感兴趣区域内。本实施例将空气-上皮分界面的邻域宽度设置为30像素，即中心线上下各15像素，参见图5所示的空气-上皮分界面ROI示意图，其中图5(a)中示意出了空气-上皮分界面ROI的上界、中心线及下界，上界和下界限定的环形区域即空气-上皮分界面ROI。

二、感兴趣区域扁平化

采用扁平化操作的目的是让后续的平均灰度差分方法获得的分界面更准确。扁平化操作以各分界面的ROI中心线和顶点纵坐标为基准，将其中的每一列像素向上移动(C_init(y)-A_y+w)像素，其中w为邻域宽度的一半，即可以得到扁平化后的ROI。参见图5(b)所示的空气-上皮分界面ROI扁平化示意图，ROI已变成矩形。需要注意的是，在采用平均灰度偏差方法获得分界面后，需要对其作逆运算，获得原图像对应的纵坐标。

三、采用平均灰度差分方法获取修正分界面

平均灰度差分计算边界公式如下：

式中j的取值范围为[1,H_ROI-1]，这里H_ROI为扁平化后ROI的高度。本发明通过计算ROI中每列的D_j最大(小)值的纵坐标来计算每列边缘的位置，对于空气-上皮分界面和上皮-鲍尔曼分界面这样的“暗到亮”型边界(从图像顶部至底部看)，可计算每列的D_j最大值的纵坐标，而对于内皮-房水分界面这样的“亮到暗”型边界，可计算每列的D_j最小值的纵坐标。该公式简单，运算量小，还可以通过构造积分图的方式来加速。考虑到噪声的影响，用该公式获取分界面后还需对其进行滤波，本实施例采用Savitzky-Golay滤波。参见图6(a)所示的采用平均灰度差分计算空气-上皮分界面并对其进行Savitzky-Golay滤波后的示意图，图6(b)为图6(a)的空气-上皮分界面进行逆扁平化后的示意图，其已对应至图6(a)中对应的位置。为了方便描述，将平均灰度差分及Savitzky-Golay滤波步骤记为步骤DoM。

在得到角膜各层的修正分界面后，还需要通过基于平抛运动的卡尔曼滤波对上述角膜各层的修正分界面进行修正和约束，从而精确地获得角膜图像各层的分界面。由于正常人眼角膜为圆形，但考虑到圆锥角膜等病眼的情况，对人眼角膜采用抛物线建模是比较合理的(相当于二阶多项式拟合)。一条抛物曲线的确定涉及到4个系数，包括顶点坐标(x_v,y_v)、曲率系数和方向，算法复杂度大，本实施例通过抛物线的内切圆来确定抛物线参数，下面首先介绍采用弦垂直平分线通过圆心的圆霍夫变换确定上述内切圆。

霍夫变换是一种重要的圆检测方法，本实施例采用的弦垂直平分线过圆心的霍夫变换，利用了圆的两个属性：不共线的三点可以确定一个圆，圆上弦垂直平分线必定通过圆心，参见图7(a)所示的垂直平分弦的线段相交于圆心的示意图，具体的霍夫变换过程如下：

设置名为KMN的三元组，该三元组包括三个点K(x_k,y_k)，M(x_M,y_M)，N(x_N,y_N)以及依次连接它们的线段KM，KN。由线段PO垂直平分线段KM并且线段QO垂直平分线段MN可知：

(y_O-y_P)×(y_K-y_M)+(x_O-x_P)×(x_K-x_M)＝0 (9)

(y_O-y_Q)×(y_M-y_N)+(x_O-x_Q)×(x_M-x_N)＝0 (10)

求解如上方程可以获得圆心：

上述两式有解的存在条件是:

(x_K-x_M)×(y_M-y_N)-(x_M-x_N)×(y_K-y_M)≠0 (13)

获取初始边界之后，可以采用上述的弦垂直平分线过圆心原理霍夫变换来获取圆C_hough。参见图7(b)所示的用于对圆心进行投票的序列等间距三元组的示意图，其中，S0至S8序列为在x方向上的等间距，令三元组KMN相邻两点的x方向距离为DIST_x像素，则K的x坐标其范围为[0，W¹-2xDIST_x-1]。霍夫变换投票分为两步，首先是圆心，然后是半径。为了更快地获取圆心和半径，本实施例从Hough空间的数字化和搜索范围上作了简化，圆心的Hough空间的投票单元其x，y方向和半径直方图其投票单元均代表1个像素，圆心的Hough空间x方向搜索范围为[Ax-100,Ax+100]，令根据先验知识得到的圆的半径约为r像素，圆心的Hough空间y方向搜索范围为[Ay+r-100,Ay+r+100]，半径直方图搜索范围为[r-100,r+100]，在圆心Hough搜索范围中求取最大值对应的x坐标和y坐标即为估计的圆心，在半径直方图搜索范围中求取最大值对应的半径即为估计的半径。上述投票可以参见图8(a)所示的霍夫变换的圆心投票直方图和图8(b)所示的霍夫变换的半径投票直方图。上述霍夫变换的结果可以参见图9所示的通过霍夫变换获得的空气-上皮分界面示意图。

从整体上看，本实施例采用的霍夫变换算法复杂度较低，原因如下：首先，由于依然存在层状结构，二维角膜分割依然可以转化为一维的定位问题；为了进一步提高效率，可按顺序等间隔地获取三元组，因此本实施例在进行霍夫变换时其数据源存储空间复杂度为一维，时间复杂度也是一维；第二，虽然两条弦垂直平分线的交点是投票在二维阵列上，但是，人眼角膜各分界面其半径在一定范围内，一旦确定角膜顶点的初始位置之后，可以获知角膜各分界面的大致圆心，这种二维累积阵列的搜索范围很小；最后，一旦确定了圆心，各边界点至圆心距离的直方图就可以确定，该直方图是一维的。

在进行基于平抛运动的卡尔曼滤波前，先证明外切于圆的抛物线，若与圆相交于切点以外的另一点，则该抛物线是唯一确定的，具体证明过程如下：

令抛物线β外切圆α于原点M(0,0)，并相交于点N(x_N,y_N),圆α的圆心为O(x_O,y_O)。可知，抛物线β和圆α的方程可写成：

(y-r)²+x²＝r² (15)

式中a为抛物线待定系数，r为圆的半径。由于抛物线和圆相交于点N，联立方程(14)、(15)并根据韦达定理可求得：

即证抛物线β满足存在性和唯一性。

因此，在对各层分界面进行霍夫变换圆拟合之后，即可以用由其内含圆限定的抛物线来对边界曲线进行卡尔曼滤波建模，以同时满足快速性和灵活性的要求。式(16)中r即霍夫变换之后圆的半径，x_N即圆与图像左右边界交点的横向坐标，结合式(14)可知，a即后续卡尔曼滤波要求的平抛运动加速度，a用来确定卡尔曼滤波模型的状态方程。

以下详述基于平抛运动的卡尔曼滤波获取角膜各层的优化分界面，参见图10(a)图10(b)所示的平抛运动示意图，可以将角膜每层分界面类比为从角膜层分界面顶点出发的两段方向相反的质点平抛运动轨迹。此时，某一质点(对应角膜层分界面上的点)以竖直方向加速度u和水平方向初速度v₀做平抛运动。此处结合角膜的形状特性，采用平抛运动原理卡尔曼滤波方程来跟踪每层上相邻的边界点。

为了逼近分界面上的点，卡尔曼滤波可以表示为如下(式17-21)，其中状态Ф和方差P能够按如下公式预测：

φ(k)＝Aφ(k-1)+Ba (17)

式中

是二维向量，y是层边界点纵坐标值，v_y是相邻边界点(质点)的速度，下标k表示当前边界点，下标(k-1)表示上一个边界点。A是状态转移矩阵，B是控制输入矩阵，^表示预测值，T表示转置。根据平抛运动方程，可知A、B分别为

式中dt表示从前一边界点至当前边界点的空间间隔，本实施例设为1，u表示质点平抛运动的竖直向下的加速度，校正方程表示如下：

式中H＝[1 0]，R是测量噪声的标准差，实际使用可设为一固定值，本实施例设为1.0，K将随着R进行调整，当R相对于

较大的时候，因子K将比较小，反之亦然。

状态Ф将会按如下公式更新：

式中Z(k)为通过步骤DoM获得的当前边界点纵坐标的值。最后，方差P将会按照如下公式更新：

总之，卡尔曼滤波通过两部分(式17-18进行预测，式19-21进行校正)来获得具有最小方差的状态。

其中，上述公式(17)，即卡尔曼滤波的状态方程φ(k)＝Aφ(k-1)+Ba可以从前面的

得到。

通过上述步骤可知，本发明实施例中分界面的获取要经过三步:首先是用平均差分方法获得初始边缘，然后是用抛物线拟合(其系数由圆来确定)，最后是卡尔曼滤波。每一步操作其精度都能够得到提高，进而可以得到角膜各层的优化分界面，参见图11所示的角膜各层的优化分界面示意图，其中自上而下依次为空气-上皮分界面，上皮-鲍尔曼分界面和内皮-房水分界面。

本发明实施例提供的从OCT角膜图像分割角膜结构的方法，首先对角膜图像进行缩放并采用横向投影和轴向投影的方式来获得顶点坐标并进一步通过平均灰度差分获取角膜各层分界面；使用弦垂直平分线过圆心原理的霍夫变换，并通过边界位置降维，分两步投票来使得霍夫变换的计算复杂度和存储复杂度降低；使用基于平抛运动原理的卡尔曼滤波来对上述各层分界面进行修正和约束，计算速度快且精度高。

实施例二

本发明实施例二提供了一种从OCT角膜图像分割角膜结构的***，参见图12所示的从OCT角膜图像分割角膜结构的***的示意图，该***包括分辨率调整模块121、顶点坐标确定模块122、初始界面确定模块123、环形感兴趣区确定模块124、矩形感兴趣区确定模块125、修正分界面确定模块126、霍夫变换圆拟合模块127和卡尔曼滤波模块128，其中，各模块的功能如下：

分辨率调整模块121，用于根据OCT角膜图像的OCT设备的分辨率，对OCT角膜图像进行缩放处理，得到横向分辨率和纵向分辨率相同的OCT等分辨率图像。顶点坐标确定模块122，用于根据OCT等分辨率图像中各像素的灰度值确定角膜顶点的横坐标和角膜各层顶点的初始纵坐标。初始界面确定模块123，用于根据角膜顶点的横坐标和角膜各层顶点的初始纵坐标，以圆拟合方式确定角膜各层的初始分界面。环形感兴趣区确定模块124，用于按照设定的邻域宽度对角膜各层的初始分界面进行边界扩充，得到角膜各层的环形感兴趣区。矩形感兴趣区确定模块125，用于对角膜各层的环形感兴趣区分别进行扁平化处理，得到角膜各层的矩形感兴趣区。修正分界面确定模块126，用于对角膜各层的矩形感兴趣区分别进行平均灰度差分处理和扁平化处理的逆处理，得到角膜各层的修正分界面。霍夫变换圆拟合模块127，用于对于角膜各层的修正界面中的每个修正分界面，基于垂直平分线过圆心原理，对每个修正分界面分别进行霍夫变换圆拟合操作，得到角膜各层的圆心和半径。卡尔曼滤波模块128，用于基于角膜各层的圆心和半径，分别对角膜各层的修正分界面进行平抛运动的卡尔曼滤波，得到角膜各层的优化分界面。

本发明实施例提供的从OCT角膜图像分割角膜结构的***，首先对角膜图像进行缩放及坐标确定，再通过平均灰度差分确定角膜各层分界面，再通过抛物线内切圆确定抛物线系数，进而采用平抛运动建模的卡尔曼滤波来对平均灰度差分方法获得的分界面进行修正和约束，可以快速并精确地获得角膜图像各层的边界。

其中，分辨率调整模块121还用于：获取OCT角膜图像的OCT设备的分辨率，其中，OCT设备的横向分辨率记为A1，轴向分辨率记为B1；获取OCT角膜图像的列数和行数，列数记为C1，行数记为D1；设置OCT等分辨率图像的列数为：C1×A1/B1，行数为D1；或者设置OCT等分辨率图像的列数为：C1，行数为D1×B1/A1；以使OCT等分辨率图像的横向分辨率和纵向分辨率相同。

其中，顶点坐标确定模块122包括：

投影单元，用于对OCT等分辨率图像分别进行轴向投影和横向投影；横坐标确定单元，用于根据轴向投影得到的灰度值大小确定角膜顶点的横坐标；中央伪影区确定单元，用于根据轴向投影得到的灰度值大小确定OCT等分辨率图像的中央伪影区；水平伪影区确定单元，用于将横向投影得到的最大灰度值对应区域，确定为OCT等分辨率图像的水平伪影区；初始纵坐标确定单元，用于基于确定出的中央伪影区和水平伪影区，获取角膜各层顶点的初始纵坐标。

进一步地，水平伪影区确定单元还用于：将OCT等分辨率图像沿横向坐标轴等宽度地分为三个区域，依次为区域I、区域Ⅱ和区域Ⅲ，根据轴向投影计算区域I和Ⅲ的平均灰度，分别记为m_I和m_III，并设定阈值T_cent为：

其中，g_max为轴向投影中的最大值，最大值对应的横坐标为P_max，ρ为预先设定的值；对轴向投影得到的灰度值，以P_max对应的位置为中心，分别向两侧逐一查找小于T_cent的灰度值，将两侧中首个查找到的灰度值对应的位置分别记为p_left和p_right；将OCT等分辨率图像中p_left和p_right之间的对应区域确定为中央伪影区。

进一步地，初始纵坐标确定单元还用于：对横向投影组成的向量进行一阶前向差分运算，将运算结果中的最大向量对应的纵坐标记为角膜的上皮层顶点的纵坐标；将OCT等分辨率图像中的中央伪影区的每列像素的灰度值设置为0，对上皮层顶点以下区域进行扁平化操作，得到OCT等分辨率图像的局部拉伸图像；对局部拉伸图像进行再次横向投影，对再次横向投影组成的向量进行一阶前向差分运算，将运算结果中的最大向量对应的纵坐标记为角膜的鲍尔曼层顶点的纵坐标；将运算结果中的最小向量对应的纵坐标记为角膜的内皮层顶点的纵坐标。

考虑到抑制水平噪声的需要，上述***还包括：

水平伪影区抑制模块，用于按照以下公式对水平伪影区进行抑制操作：

其中，f(x,y)为OCT等分辨率图像中像素(x,y)处的灰度值，W为OCT等分辨率图像的横向长度、H为OCT等分辨率图像的纵向长度，x与y为非负整数，取值范围分别为[0,W-1]与[0,H-1]其中min与max分别为

的最小值与最大值。

进一步地，修正分界面确定模块126还用于：采用以下公式分别计算角膜各层的矩形感兴趣区的界面：

其中，j的取值范围为[1,H_ROI-1]，H_ROI为扁平化后矩形感兴趣区的高度；f(x,y)为矩形感兴趣区中像素(x,y)处的灰度值；角膜各层的矩形感兴趣区包括：空气-上皮分界面对应的矩形感兴趣区、上皮层-鲍尔曼层分界面对应的矩形感兴趣区和内皮层-房水分界面对应的矩形感兴趣区；将空气-上皮分界面对应的矩形感兴趣区计算出的每列D_j中的最大值作为上皮层界面在该列的纵坐标；将上皮层-鲍尔曼层分界面对应的矩形感兴趣区计算出的每列D_j中的最大值作为鲍尔曼层界面在该列的纵坐标；将内皮层-房水分界面对应的矩形感兴趣区计算出的每列D_j中的最小值作为内皮层界面在该列的纵坐标。

进一步地，霍夫变换圆拟合模块127还用于：对角膜各层的每个修正分界面，均进行如下操作得到角膜各层的圆心和半径：

按照设定的第一间隔，依次从当前修正分界面上选取三个像素点K、M和N，连接KM和MN；将KM和MN各自垂直平分线的交点确定为当前圆心；按顺序选取K、M和N，以及圆心的确定过程，直至当前修正分界面上各像素点均被选取过为止，将重合次数最多的圆心作为当前修正分界面的最终圆心；以最终圆心为起点，按照设定的第二间隔，从当前修正分界面上逐一选取像素点，计算选取的像素点到最终圆心的距离，将距离作为当前半径；将最多重复次数的半径确定为当前修正分界面的半径。

进一步地，卡尔曼滤波模块158还用于：分别以角膜各层的顶点为原点建立坐标系，按照下述公式获取角膜各层的优化分界面：

y＝ax²/2；其中，

x_N为各个修正分界面对应的圆与OCT等分辨率图像左边界或右边界交点的横向坐标，r为修正分界面对应的圆的半径；

按照如下方式设置卡尔曼滤波模型中的状态参数Ф和方差P：

φ(k)＝Aφ(k-1)+Ba；

其中，

是二维向量，y是当前角膜层的优化分界面上各点的纵坐标值，v_y是相邻边界点的速度，下标k表示当前边界点，下标(k-1)表示上一个边界点；A是状态转移矩阵，B是控制输入矩阵，–表示预测值，T表示转置；

式中dt表示从前一边界点至当前边界点的空间间隔；

构建卡尔曼滤波模型中的校正方程如下：

其中，H＝[10]，R是预先设定的测量噪声的标准差，K为随着R进行调整的参数；

设置状态参数Ф按如下公式进行更新：

其中，Z_k为当前角膜层的修正分界面上的k点的纵坐标值；

设置方差P按照如下公式更新：

将更新后的Ф和P代入上述校正方程，直至方差P趋近于最小方差值。

本发明实施例提供的从OCT角膜图像分割角膜结构的***，对角膜图像进行缩放并采用横向投影和轴向投影的方式来获得顶点坐标并进一步通过平均灰度差分获取角膜各层分界面；使用弦垂直平分线过圆心原理的霍夫变换，并通过边界位置降维，分两步投票来使得霍夫变换的计算复杂度和存储复杂度降低；使用基于平抛运动原理的卡尔曼滤波来对上述各层分界面进行修正和约束，计算速度快且精度高。

本发明实施例所提供的***，其实现原理及产生的技术效果和前述方法实施例相同，为简要描述，***实施例部分未提及之处，可参考前述方法实施例中相应内容。

以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应所述以权利要求的保护范围为准。

Claims (10)

1.一种从OCT角膜图像分割角膜结构的方法，其特征在于，包括：

根据OCT角膜图像的OCT设备的分辨率，对OCT角膜图像进行缩放处理，得到横向分辨率和纵向分辨率相同的OCT等分辨率图像；其中，所述OCT角膜图像为所述OCT设备生成的二维灰度图；

根据所述OCT等分辨率图像中各像素的灰度值确定角膜顶点的横坐标和角膜各层顶点的初始纵坐标，其中，所述角膜各层顶点包括：上皮层顶点、鲍尔曼层顶点和内皮层顶点；

根据所述角膜顶点的横坐标和所述角膜各层顶点的初始纵坐标，以圆拟合方式确定所述角膜各层的初始分界面；

按照设定的邻域宽度对所述角膜各层的初始分界面进行边界扩充，得到所述角膜各层的环形感兴趣区；

对所述角膜各层的环形感兴趣区分别进行扁平化处理，得到所述角膜各层的矩形感兴趣区；

对所述角膜各层的矩形感兴趣区分别进行平均灰度差分处理和所述扁平化处理的逆处理，得到所述角膜各层的修正分界面；

对于所述角膜各层的修正界面中的每个修正分界面，基于垂直平分线过圆心原理，对每个所述修正分界面分别进行霍夫变换圆拟合操作，得到所述角膜各层的圆心和半径；

基于所述角膜各层的圆心和半径，分别对所述角膜各层的修正分界面进行平抛运动的卡尔曼滤波，得到所述角膜各层的优化分界面。

2.根据权利要求1所述的方法，其特征在于，根据OCT角膜图像的OCT设备的分辨率，对OCT角膜图像进行缩放处理，得到横向分辨率和纵向分辨率相同的OCT等分辨率图像包括：

获取OCT角膜图像的OCT设备的分辨率，其中，所述OCT设备的横向分辨率记为A1，轴向分辨率记为B1；

获取所述OCT角膜图像的列数和行数，所述列数记为C1，所述行数记为D1；

设置OCT等分辨率图像的列数为：C1×A1/B1，行数为D1；或者设置OCT等分辨率图像的列数为：C1，行数为D1×B1/A1；以使所述OCT等分辨率图像的横向分辨率和纵向分辨率相同。

3.根据权利要求1所述的方法，其特征在于，根据所述OCT等分辨率图像中各像素的灰度值确定角膜顶点的横坐标和角膜各层顶点的初始纵坐标包括：

对所述OCT等分辨率图像分别进行轴向投影和横向投影；

根据所述轴向投影得到的灰度值大小确定角膜顶点的横坐标；

根据所述轴向投影得到的灰度值大小确定所述OCT等分辨率图像的中央伪影区；

将所述横向投影得到的最大灰度值对应区域，确定为所述OCT等分辨率图像的水平伪影区；

基于确定出的中央伪影区和水平伪影区，获取角膜各层顶点的初始纵坐标。

4.根据权利要求3所述的方法，其特征在于，根据所述轴向投影得到的灰度值大小确定所述OCT等分辨率图像的中央伪影区包括：

将所述OCT等分辨率图像沿横向坐标轴等宽度地分为三个区域，依次为区域I、区域Ⅱ和区域Ⅲ，根据所述轴向投影计算区域I和Ⅲ的平均灰度，分别记为m_I和m_III，并设定阈值T_cent为：

其中，g_max为所述轴向投影中的最大值，所述最大值对应的横坐标为P_max，ρ为预先设定的值；

对所述轴向投影得到的灰度值，以P_max对应的位置为中心，分别向两侧逐一查找小于T_cent的灰度值，将两侧中首个查找到的灰度值对应的位置分别记为p_left和p_right；

将所述OCT等分辨率图像中p_left和p_right之间的对应区域确定为中央伪影区。

5.根据权利要求3所述的方法，其特征在于，基于确定出的中央伪影区和水平伪影区，获取角膜各层顶点的初始纵坐标包括：

对所述横向投影组成的向量进行一阶前向差分运算，将运算结果中的最大向量对应的纵坐标记为所述角膜的上皮层顶点的纵坐标；

将所述OCT等分辨率图像中的中央伪影区的每列像素的灰度值设置为0，对所述上皮层顶点以下区域进行扁平化操作，得到所述OCT等分辨率图像的局部拉伸图像；

对所述局部拉伸图像进行再次横向投影，对所述再次横向投影组成的向量进行一阶前向差分运算，将运算结果中的最大向量对应的纵坐标记为所述角膜的鲍尔曼层顶点的纵坐标；将运算结果中的最小向量对应的纵坐标记为所述角膜的内皮层顶点的纵坐标。

6.根据权利要求3所述的方法，其特征在于，确定出所述OCT等分辨率图像的水平伪影区之后，所述方法还包括：

按照以下公式对所述水平伪影区进行抑制操作：

其中，f(x,y)为所述OCT等分辨率图像中像素(x,y)处的灰度值，W为所述OCT等分辨率图像的横向长度、H为所述OCT等分辨率图像的纵向长度，x与y为非负整数，取值范围分别为[0,W-1]与[0,H-1]其中min与max分别为

的最小值与最大值。

7.根据权利要求1所述的方法，其特征在于，对所述角膜各层的矩形感兴趣区分别进行平均灰度差分处理包括：

采用以下公式分别计算所述角膜各层的矩形感兴趣区的界面：

其中，j的取值范围为[1,H_ROI-1]，H_ROI为扁平化后矩形感兴趣区的高度；f(x,y)为所述矩形感兴趣区中像素(x,y)处的灰度值；所述角膜各层的矩形感兴趣区包括：空气-上皮分界面对应的矩形感兴趣区、上皮层-鲍尔曼层分界面对应的矩形感兴趣区和内皮层-房水分界面对应的矩形感兴趣区；

将空气-上皮分界面对应的矩形感兴趣区计算出的每列D_j中的最大值作为上皮层界面在该列的纵坐标；

将上皮层-鲍尔曼层分界面对应的矩形感兴趣区计算出的每列D_j中的最大值作为鲍尔曼层界面在该列的纵坐标；

将内皮层-房水分界面对应的矩形感兴趣区计算出的每列D_j中的最小值作为内皮层界面在该列的纵坐标。

8.根据权利要求1所述的方法，其特征在于，基于垂直平分线过圆心原理，对每个所述修正分界面分别进行霍夫变换圆拟合操作，得到所述角膜各层的圆心和半径包括：

对所述角膜各层的每个修正分界面，均进行如下操作得到所述角膜各层的圆心和半径：

按照设定的第一间隔，依次从当前修正分界面上选取三个像素点K、M和N，连接KM和MN；将KM和MN各自垂直平分线的交点确定为当前圆心；

按顺序选取K、M和N，以及圆心的确定过程，直至所述当前修正分界面上各像素点均被选取过为止，将重合次数最多的圆心作为所述当前修正分界面的最终圆心；

以所述最终圆心为起点，按照设定的第二间隔，从所述当前修正分界面上逐一选取像素点，计算选取的像素点到所述最终圆心的距离，将所述距离作为当前半径；将最多重复次数的半径确定为所述当前修正分界面的半径。

9.根据权利要求1所述的方法，其特征在于，基于所述角膜各层的圆心和半径，分别对所述角膜各层的修正分界面进行平抛运动的卡尔曼滤波，得到所述角膜各层的优化分界面，包括：

分别以所述角膜各层的顶点为原点建立坐标系，设置下述抛物线方程以获取所述角膜各层的优化分界面：

y＝ax²/2,

其中，

a为抛物线方程的系数，该系数用来确定卡尔曼滤波模型的状态方程，x_N为各个所述修正分界面对应的圆与所述OCT等分辨率图像左边界或右边界交点的横向坐标，r为所述修正分界面对应的圆的半径；

按照如下方式设置卡尔曼滤波模型中的状态参数Ф和方差P：

φ(k)＝Aφ(k-1)+Ba

其中，

是二维向量，y是当前角膜层的优化分界面上各点的纵坐标值，v_y是相邻边界点的速度，下标k表示当前边界点，下标(k-1)表示上一个边界点；A是状态转移矩阵，B是控制输入矩阵，^表示预测值，T表示转置；

式中dt表示从前一边界点至当前边界点的空间间隔；

构建卡尔曼滤波模型中的校正方程如下：

其中，H＝[1 0]，R是预先设定的测量噪声的标准差，K将随着R和

进行更新；

设置状态参数Ф按如下公式进行更新：

其中，Z(k)为所述当前角膜层的修正分界面上的当前点的纵坐标值；

设置方差P按照如下公式更新：

将更新后的Ф和P代入上述校正方程，直至方差P趋近于最小方差值。

10.一种从OCT角膜图像分割角膜结构的***，其特征在于，包括：

分辨率调整模块，用于根据OCT角膜图像的OCT设备的分辨率，对OCT角膜图像进行缩放处理，得到横向分辨率和纵向分辨率相同的OCT等分辨率图像；其中，所述OCT角膜图像为所述OCT设备生成的二维灰度图；

顶点坐标确定模块，用于根据所述OCT等分辨率图像中各像素的灰度值确定角膜顶点的横坐标和角膜各层顶点的初始纵坐标，其中，所述角膜各层顶点包括：上皮层顶点、鲍尔曼层顶点和内皮层顶点；

初始界面确定模块，用于根据所述角膜顶点的横坐标和所述角膜各层顶点的初始纵坐标，以圆拟合方式确定所述角膜各层的初始分界面；

环形感兴趣区确定模块，用于按照设定的邻域宽度对所述角膜各层的初始分界面进行边界扩充，得到所述角膜各层的环形感兴趣区；

矩形感兴趣区确定模块，用于对所述角膜各层的环形感兴趣区分别进行扁平化处理，得到所述角膜各层的矩形感兴趣区；

修正分界面确定模块，用于对所述角膜各层的矩形感兴趣区分别进行平均灰度差分处理和所述扁平化处理的逆处理，得到所述角膜各层的修正分界面；

霍夫变换圆拟合模块，用于对于所述角膜各层的修正界面中的每个修正分界面，基于垂直平分线过圆心原理，对每个所述修正分界面分别进行霍夫变换圆拟合操作，得到所述角膜各层的圆心和半径；

卡尔曼滤波模块，用于基于所述角膜各层的圆心和半径，分别对所述角膜各层的修正分界面进行平抛运动的卡尔曼滤波，得到所述角膜各层的优化分界面。

Images