基于回溯最短路径算法的相干光断层图像的层分割方法

摘要

本发明公开了一种基于回溯最短路径算法的相干光断层图像的层分割方法，包括步骤：S1、对原始图像进行去噪处理；S2、对去除噪声后的图像，计算图中所有节点权重，形成权重图；S3、采用自动初始化端点的方法在权重图上定义最短路径算法的起点和终点；S4、在权重图上用带回溯的最短路径算法计算并寻找起点与终点间的最小累计代价路径；S5、在得到的最小累计代价路径中去除自定义端点时添加的辅助列，得到相干光断层图像中其中一层的层分割结果；限制图像搜索区域，重复上述步骤S3‑S5，最终得到相干光断层图像中所有层的分割结果。该分割方法避免了传统方法中的短路问题且适用于正常与受病变影响相干断层图像的分割，具有高准确性和灵活性。

权利要求

1.一种基于回溯的最短路径算法的相干光断层图像的层分割方法，其特征在于，包括如下步骤：S1、对原始图像进行去噪处理；

S2、对去除噪声后的图像，计算图中所有节点权重，形成权重图；

S3、采用自动初始化端点的方法在权重图上定义最短路径算法的起点和终点；

S4、在权重图上用带回溯的最短路径算法计算并寻找起点与终点间的最小累计代价路径；在权重图上用带回溯的最短路径算法计算并寻找起点与终点间的最小累计代价路径采用如下公式：其中，N表示当前的图像节点， 表示路径上从起始节点start到当前图像节点N的节点权重值累加和， 表示路径上从起始节点start到图像节点N-l的节点权重值累加和， 表示当前图像节点向起始节点方向回溯长度l到前一个点N-l+1的累计代价；

节点i的代价Wi的计算公式如下：

Wi＝r*Wigray+s*Wigrad+t*Wigabor*Dif其中，Wigray为节点的图像灰度值，Wigrad为梯度权值，Wigabor为节点的Gabor变换后的最大响应值权重，r，s，t是三个常数参数，r+s+t＝1，Dif表示路径的方向偏差；

S5、在得到的最小累计代价路径中去除自定义端点时添加的辅助列，得到相干光断层图像中其中一层的层分割结果；

限制图像搜索区域，重复上述步骤S3-S5，最终得到相干光断层图像中所有层的分割结果。

2.根据权利要求1所述的方法，其特征在于，步骤S1中，采用BM3D算法进行去噪处理。

3.根据权利要求1所述的方法，其特征在于，步骤S2中，所述节点权重计算包括灰度权重、梯度权重和图像Gabor变换后的最大响应值权重。

说明书

基于回溯最短路径算法的相干光断层图像的层分割方法

技术领域

[0001] 本发明涉及相干光断层图像处理技术领域，特别是涉及一种基于回溯最短路径算法的相干光断层图像的层分割方法。

背景技术

[0002] OCT(Optical Coherence Tomography)，又称光学相干断层扫描成像，是一种新的层析成像技术，具有高分辨率，在生物内部结构上断层成像的特点。OCT之所以能够获得近显微分辨率的次表面图像，是因为它采用光波作为探测手段，而不是声波。光波的波长比声波的波长短，传播速度比声波快。OCT是基于低相干干涉技术发展的，近几十年来取得了巨大的进步。一代OCT又被称为时域OCT(TD-OCT)，每秒采集400张A扫描图，轴向分辨率为8微米到10微米。2011年出现了谱域OCT(SD-OCT)。目前商用的谱域OCT设备成像速度比时域OCT设备快10到100倍，轴向分辨率范围为3到7微米。最近又出现了扫源OCT(SS-OCT)，每秒钟可以采集超过十万张的A扫描图，轴向分辨率可以达到3微米。

[0003] OCT可以对生物组织实时地成像，而不需要制备样品。它不涉及到离子辐射，对人体健康没有危害。OCT在许多领域得到迅速的发展，通过使用低相干干涉法，OCT可以从生物样本的反向散射剖面中提供高分辨率的横截面图像。这使得在疾病进展期间或在患者认识到视力变化(视敏度)之前，可以观察到视网膜的病理和形态学变化。在过去的二十年中，OCT已经成为一种成熟的影像学方法，被眼科医生广泛应用于视网膜和光学神经疾病的诊断。视网膜层的畸变是专家诊断年龄相关性黄斑变性，黄斑水肿和黄斑裂孔的重要信号。然而，手工层分割通常是一个耗时且可重复度差的过程。因此，OCT图像的自动或半自动分割方法是必要的。

[0004] 目前出现了一些自动或半自动的OCT层分割方法，这些方法大致可分为三类:基于A扫描的方法、基于B扫描的方法和基于三维体数据的方法。基于A扫描的方法可以在每一个A扫描剖面上检测到边界上的强度峰值点或谷值点，然后用模型拟合技术将检测到的点连接起来，形成一个平滑和连续的边界。然而，基于A扫描的方法大多不考虑探测点周围的信息，容易受到噪声的影响，因此，基于A扫描的方法的准确性和鲁棒性是有限的。基于B扫描的常用方法包括活动轮廓法、基于最短路径的图搜索方法以及统计学模型的方法。B扫描方法在一般情况下优于A扫描方法，然而，它们仍有可能无法检测到病变的视网膜结构。现有的基于体数据的分割方法主要采用基于3D图的方法和模式识别方法，但这些方法一般计算成本较高，模式识别方法通常需要包括由专家手动分割的训练数据来学习一种可行的分类模型。这些模式识别方法也在准确性和效率方面受到影响。尽管适用于正常图像或有限病理的图像算法已经取得了很大的进展，但对具有明显病理变化的视网膜图像层分割仍然是一个挑战。例如，在黄斑裂孔中，当中心区域出现严重变形时，上述方法将不能准确地进行分割。

[0005] 用于求解图中最短路径问题的算法称为“最短路径算法”，最常用的最短路径算法有:Dijkstra算法、Bellman-Ford算法和Floyd-Warshall算法。在最短路径法中，它通过在两个给定的点(起点和终点)之间的连通路径，以最小的累计代价来提取图像中的曲线状结构。但是传统最短路径算法具有短路问题，本发明提出改进的带回溯最短路径算法，可以有效避免短路径问题，适用于病变视网膜层结构分割图像。

发明内容

[0006] 本发明为了克服以上技术的不足，提供了一种基于回溯最短路径算法的相干光断层图像的层分割方法。

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

[0008] 一种基于回溯最短路径算法的相干光断层图像的层分割方法，包括如下步骤：

[0009] S1、对原始图像进行去噪处理；

[0010] S2、对去除噪声后的图像，计算图中所有节点权重，形成权重图；

[0011] S3、采用自动初始化端点的方法在权重图上定义最短路径算法的起点和终点；

[0012] S4、在权重图上用带回溯的最短路径算法计算并寻找起点与终点间的最小累计代价路径；

[0013] S5、在得到的最小累计代价路径中去除自定义端点时添加的辅助列，得到相干光断层图像中其中一层的层分割结果；

[0014] 限制图像搜索区域，重复上述步骤S3-S5，最终得到相干光断层图像中所有层的分割结果。

[0015] 本发明优选的，步骤S1中，采用BM3D算法进行去噪处理，该方法快捷方便。

[0016] 本发明优选的，步骤S2中，所述节点权重计算包括灰度权重、梯度权重和图像Gabor变换后的最大响应值权重。

[0017] 本发明优选的，步骤S4中，在权重图上用带回溯的最短路径算法计算并寻找起点与终点间的最小累计代价路径采用如下公式：

[0018]

[0019] 其中，N表示当前的图像节点， 表示当前图像节点向起始节点方向回溯长度l到前一个点N-l+1的累计代价；；

[0020] 节点Wi的代价计算公式如下：

[0021] Wi＝r*Wigray+s*Wigrad+t*Wigabor*Dif

[0022] 其中，Wigray为节点的图像灰度值，Wigrad为梯度权值，Wigabor为节点的Gabor最大响应值，r，s，t是三个常数参数，r+s+t＝1。

[0023] 本发明的有益效果是：

[0024] 本发明的方法首先将病变的OCT图像去除噪声，基于降噪后的图像，计算图中所有节点权重图，然后通过自动初始化端点方法定义起点和终点，最后在权重图上用带回溯的最短路径算法计算起点和终点间的最小累计代价路径作为视网膜层分割结果。本发明的分割方法具有如下优点：

[0025] (1)高准确性。本方法使用带回溯的最短路径算法，首次在该技术领域使用，进行了病变相干断层图像层结构分割，避免了传统方法中的短路问题，与现有的方法相比，得到更好的分割效果和性能。

[0026] (2)灵活性。本方法同样适用于正常相干断层图像的层结构分割，可以在正常图像上得到较好的分割结果图，具有较高的灵活性和适应性。

附图说明

[0027] 图1为本发明实施例的流程示意图。

[0028] 图2为视网膜层结构示意图。

[0029] 图3(a)为Gabor变换前原图，图3(b)、图3(c)、图3(d)分别是沿30度方向、水平方向、垂直方向的响应图，图3(e)是Gabor变换后的最大响应图，图3(f)是Gabor后的方向信息图。

[0030] 图4是分割前原图。

[0031] 图5是图4的层分割结果图。

[0032] 图6是自动初始化起点终点的示意图。

具体实施方式

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

[0034] 本发明公开了一种基于回溯最短路径算法的相干光断层图像的层分割方法，下面以病变视网膜层为例并结合附图1-6说明其具体实施过程。

[0035] 步骤一、对原始图像进行去噪处理

[0036] 本发明所用到的视网膜图像数据带有噪声，会影响最终的图像分割结果，现有很多图像去噪方法，本实施例采用BM3D算法进行快捷方便的去噪处理，得到去噪后的图像。

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

[0038] 步骤二、计算节点权重得到权重图

[0039] 经步骤一得到去噪后的图像后，第二步需计算去噪后图像上所有节点权重，在本实施例中，采用图像强度、梯度信息和Gabor方向信息进行节点权重计算，具体如下：

[0040] 首先，图像的强度权重计算。图像强度表示单通道图像像素的强度，在我们的数据集中，图像是灰度图像，所以强度是像素的灰度，在视网膜OCT图像中，本实施例强度权重的计算公式如下：

[0041]

[0042] 其中，I(i,j)表示原始图像中节点的图像灰度值，i,j为图像中像素的位置，k是一个常数值，设为1。

[0043] 其次，梯度信息的权重计算。梯度图像是由原始图像通过与滤波器卷积得到的，图像梯度可以用来提取图像的边缘信息，在本实施例中，视网膜OCT图像中的层可以作为边缘信息，因此要用到梯度信息，图像梯度权重的计算公式如下：

[0044]

[0045] dx(i,j)＝I(i+1,j)-I(i,j)  (3)

[0046]

[0047] 其中，I是图像，i,j是像素的坐标，考虑到边缘有两种：暗到亮的边缘和亮到暗的边缘，我们计算dy(i,j)如上所示。

[0048] 最后，Gabor方向信息的权重计算。Gabor变换是短时傅里叶变换的一个特例，当信号随时间变化时，通常用来确定信号的正弦频率和相位内容，Gabor变换可用于提取不同方向的特征。在本方法中，Gabor变换应用于每个像素，以计算其在不同方向上的响应值，然后使用Gabor最大响应值作为权重进行计算。

[0049] 当计算一个节点的权重时，结合上面介绍的三种权重计算方法如下：

[0050] Wi＝r*Wigray+s*Wigrad+t*Wigabor*Dif  (5)

[0051] 其中，Wigray为节点的图像灰度权重，Wigrad为梯度权值，Wigabor为节点的最大响应值，r，s，t是三个常数参数，它们的和r+s+t＝1。

[0052] 步骤三、自动初始化算法传播起点和终点

[0053] 图6显示了一个带有自动初始化技术的示例图像。首先，本方法在步骤二得到的权重图左、右两侧各添加一列，这两列的权重值设为整个权重图的最小值；然后，起点和终点自动初始化为左上角节点和右下角节点，如图6所示，最左一列和最右一列为添加的辅助权重最小值列。

[0054] 步骤四、在权重图上用带回溯的最短路径算法计算最小累计代价路径

[0055] 由于存在短路问题，原有的最短路径算法不能准确地分割出畸变后的视网膜层结构，因此，本实施例采用回溯法来解决这个问题。在权重图起点与终点间计算累计代价时，从当前节点向起点方向回溯一定长度，以回溯的短边累计成本代替节点代价。这样做，可以避免短路的问题，公式如下：

[0056]

[0057] 其中，N是当前的图像节点， 表示当前图像节点回溯长度l到前一个点的累计代价。

[0058] 此外，本实施方法还利用Gabor方向信息来定义一个新的权重分量来限制算法的传播方向。在OCT图像的层分割任务中，层的方向通常是水平的，如图2所示，而这一先验信息可以被利用来提高分割方法的鲁棒性。如图3(f)所示，当gabor变换应用于图像时，每个节点都有一个主方向值，本方法中用这个方向值来限制传播方向。在算法传播过程中，我们计算了回溯后短边的平均方向值，并将这个平均方向作为传播趋势，将当前节点的方向值将与此平均方向值进行比较，确定下一个节点的传播。

[0059] 步骤五、在得到的最小累计代价路径中去除自定义端点时添加的辅助列，得到相干光断层图像中其中一层的层分割结果

[0060] 如图6所示，在添加辅助列后的权重图上，线条表示使用本实施例方法后的示例结果，然后将添加的左右两列辅助列删除，得到病变视网膜层结构中ILM层的分割结果。

[0061] 分割出ILM层之后，限制图像搜索区域，重复上述步骤三到五，得到如图2所示其余各层的分割结果。

[0062] 本发明的描述中，术语“水平”、“垂直”、“上”、“下”、“左”、“右”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了描述本发明和简化描述，而不是指示或暗示所指的结构必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。

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