基于曲线群匹配的OCT眼底图像半自动分割方法和装置

技术领域

本发明涉及医学图像处理领域，特别是一种基于曲线群匹配的OCT眼底图像半自动分割方法和装置。

背景技术

光学相干层析成像(Optical Coherence Tomography,OCT)技术是第一种基于低相干光干涉原理的生物组织微观医学成像新技术，其理论基础是早期的白光干涉测量法。因其高分辨率、无损伤和非入侵性的特点，在临床医学领域具有巨大的应用价值。通过OCT检查获得视网膜图像可以清晰展示眼底结构与功能信息。大量研究表明，视网膜神经纤维层(retinal never fiber layer,RNFL)厚度的变化可作为青光眼、黄斑病变等神经变性疾病发病的重要指标。OCT可以获取视网膜神经纤维层在各处的厚度，并能够检测出传统视网膜厚度分析仪所不能识别出的弥漫性萎缩，OCT提供RNFL厚度测量，能够及早发现视网膜神经纤维的弥漫性变薄，减少青光眼的发病率。通过对OCT视网膜图像的分析处理获得视网膜厚度的数据和地形图，可以测量出黄斑区周围的视网膜后度，有助于深入了解视网膜黄斑异常状态，监测黄斑裂孔、黄斑水肿、老年性黄斑病变进展。OCT提供的视网膜层状结构图像提高了医学对视网膜的结构和功能的理解，在揭示视网膜疾病的发病机理，确定新型治疗方案，监测疾病治疗效果等方面，发挥越来越重要的作用，已经成为眼科检查的常规手段。

近几年，计算机辅助分割吸引了大量的研究注意力，因为它不仅能够提供可重复的，定量的和客观的图像分割结果，而且减少医生进行OCT视网膜图像分割时所需的时间和精力。因此科研人员研究了大量的分割算法应用于OCT眼底图像的分割。其中，大多数算法是通过探索视网膜每层区域内以及层与层边界处的灰度和几何特性进行视网膜层分割，比如：灰度变化测量、像素分类、马尔科夫边界建模、边缘检测、纹理和形状分析、机器学习技术和图像配准等技术。在这些算法中，基于图论的分割方法已经被证明是成功的，并且有研究指出它在各种医学图像分割的应用中是优于经典的最短路径算法的。但是这类方法在处理层间对比度低、有严重视网膜畸形的真实OCT视网膜图像时产生的结果大多数时候精度是有限的。这种情况下，部分研究者通过在其中加入人工交互部分对原有算法进行改进，而人为参与分割的过程不仅耗时还可能导致结果不客观。此外，大部分现有的基于计算机的视网膜层分割算法在不进行技术重新设计或实验重新训练的情况下，不能分割不同的视网膜层或者来自不同物种的视网膜层，例如人类与犬类的视网膜。这种局限性将相关技术的应用限制在一个非常特定的领域。

如何提高OST眼底图像分割的效率，以及增强分割算法的适应性，是本领域技术人员目前迫切解决的技术问题。

发明内容

为克服上述现有技术的不足，本发明结合了动态时间规划和多目标特征点匹配技术，提出了一种新的基于曲线群匹配的半自动化分割模型，该模型在不需要进行任何额外工作的情况下能够对眼底OCT图像中的任意视网膜层以及部分病变导致的层状结构都能够进行准确分割，不受大部分图像中中央凹处各层结构区别不明显这一现象的影响。具体来说，我们将OCT图像每一列灰度值由上到下的强度变化生成一条曲线，并且将来自不同列的曲线以分组方式进行两两匹配。当曲线匹配生成成对曲线之间的空间对应关系后，我们手动指定任何视网膜层上的任一边界点便能够确定视网膜层的边界。

为实现上述目的，本发明采用如下技术方案：

一种基于曲线群匹配的OCT眼底图像半自动分割方法，包括以下步骤：

(1)接收一幅m行n列的眼底图像并进行预处理；

(2)确定所述眼底图像每一列的灰度值向量曲线f_α＝[x₁,x₂,x₃,…,x_m]^T，α∈{1,n}，x_i∈{0,255}，计算灰度值向量曲线f_α中每个像素点的描述子，得到所述眼底图像的描述子矩阵；

(3)基于所述描述子矩阵，通过动态时间规划算法将所有灰度值向量曲线进行两两匹配，得到成对曲线像素点之间的空间对应关系矩阵；

(4)接收用户手动指定的像素点，根据所述点的坐标和所述空间对应关系矩阵，找到与所述点对应的所有像素点，拟合得到平滑的分割线。

进一步地，所述步骤(1)的预处理包括图像滤波和对比度增强。

进一步地，所述步骤(2)像素点的描述子的计算方法为：

对于灰度值向量曲线f_α＝[x₁,x₂,x₃,…,x_m]^T，α∈{1,n}，x_i∈{0,255}，针对任一元素x_i，取距离该点最近的l个点构成子序列p_i，将p_i分成s个可重叠的小区间，每个小区间的维度小于l；

用最小二乘法在每个小区间内拟合直线，将子序列p_i内所有小区间的直线斜率整合到一起，得到像素点x_i的s维的向量描述子d_i＝f_slope(p_i)，表示x_i的局部曲线特征。

进一步地，所述步骤(3)包括：

采用欧氏距离计算任意两条曲线上各向量描述子的相似度，得到相似度矩阵；

基于动态时间规划模型求解两条曲线上点的索引对应关系矩阵P；

所述动态时间规划函数为：

其中，P＝[p^α,p^β]^T∈R^2×k表示两个描述子序列索引的对应关系，k表示对齐两个各描述子所需要的步数，当曲线D_α中的第i个描述子与曲线D_β中的第j个描述子在第t步对应时，则所述模型的求解是通过路径规划策略寻找相似度矩阵中相似度最小的路径；

根据索引对应关系矩阵P，计算得到对应矩阵L_αβ∈R^m×m，当且仅当存在时，L_αβ中第(i,j)元素为1；把每两条曲线的对应矩阵计算出来后存储到一个连接矩阵L∈{0,1}^mn×mn中，得到成对曲线之间的空间对应关系矩阵，其中，当α＝β时，L_αβ＝1。

进一步地，所述路径规划策略为：

其中，L^*(P_r)表示到第r步所需要的损失函数，ω表示连续步骤间的转化方式，约束条件包括：

(1)边界条件：所有规划路径都是从点(1,1)开始到(m,m)结束，即：从步骤P₁＝[1,1]^T开始到步骤P_k＝[m,m]^T结束；

(2)单一性：如果r₁≤r₂，则p_r1-p_r2≤0；

(3)连续性：P_r-P_r-1∈{[0,1]^T,[1,0]^T,[1,1]^T}。

进一步地，所述步骤(3)得到空间对应关系矩阵后还包括：对所述空间对应关系矩阵进行像素点的匹配优化。

进一步地，所述匹配优化包括：

构建优化目标函数：

其中，<·，·>表示内积，X_αβ表示L_αβ优化后的矩阵，α是平衡稀疏结构程度的正则项参数，λ为核函数的权重，核范数||X||_*最小化X的秩，X需满足以下约束：(1)自身匹配被记为单位矩阵，(2)X必须是对称矩阵，(3)X中的值必须在[0,1]之间；

用Z表示α1-L，令X＝AB^T，u＜N，A,B∈R^N×u，优化函数变形为：

其中表示满足上述三个约束条件的所有矩阵的集合；

采用交替方向乘子算法解决优化函数最小化问题。

进一步地，所述步骤(4)采用三次B样条插值进行分割线拟合。

根据本发明的第二目的，本发明还提供了一种基于曲线群匹配的OCT眼底图像半自动分割装置，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，所述处理器执行所述程序时实现所述的基于曲线群匹配的OCT眼底图像半自动分割方法。

根据本发明的第三目的，本发明还提供了一种计算机可读存储介质，其上存储有计算机程序，该程序被处理器执行时执行所述的基于曲线群匹配的OCT眼底图像半自动分割方法。

本发明的有益效果

1、在处理效果上，本发明首次提出基于曲线匹配的OCT眼底图像分割算法，并结合基于全局的优化算法。该模型能够从视网膜OCT图像中可靠精确的提取视神经纤维层轮廓并对其进行高效分割的同时，在分割对比度较低纹理结构不明显的中央凹结构是也表现出良好的性能。

2、在适用性和扩展性上，曲线匹配框架能够适用于大部分包含层状结构的图形分割情形，本方法不需要定义、提取整体特征，对局部细节的分割表现了良好的性能，能够进行健康视网膜图像以及大多数视网膜病变图像的分割；此外该模型可以进一步扩展维度进行3D视网膜OCT图像的分割。

附图说明

构成本申请的一部分的说明书附图用来提供对本申请的进一步理解，本申请的示意性实施例及其说明用于解释本申请，并不构成对本申请的不当限定。

图1是本发明基于曲线群匹配的视网膜OCT图像分割模型的流程图；

图2是本模型对健康视网膜OCT图像进行预处理的结果图；

图3是本模型对有疾病的视网膜OCT图像进行预处理的结果图；

图4是利用动态时间规划模型匹配两个灰度值曲线的一个实例；

图5是动态时间规划模型进行曲线匹配时规划路径示意图；

图6是本发明在色素上皮层脱离的视网膜图像上分割结果的示意图；

图7是本发明在神经纤维层下积液的视网膜OCT图像上分割结果的示意图。

具体实施方式

应该指出，以下详细说明都是例示性的，旨在对本申请提供进一步的说明。除非另有指明，本文使用的所有技术和科学术语具有与本申请所属技术领域的普通技术人员通常理解的相同含义。

需要注意的是，这里所使用的术语仅是为了描述具体实施方式，而非意图限制根据本申请的示例性实施方式。如在这里所使用的，除非上下文另外明确指出，否则单数形式也意图包括复数形式，此外，还应当理解的是，当在本说明书中使用术语“包含”和/或“包括”时，其指明存在特征、步骤、操作、器件、组件和/或它们的组合。

在不冲突的情况下，本申请中的实施例及实施例中的特征可以相互组合。

本发明可以总结为曲线两两匹配和根据匹配结果进行整体优化两部分：首先利用动态时间规划函数对图像中每一列提取的灰度值曲线进行两两匹配，得到所有曲线两两匹配的结果，再利用循环一致性将该匹配结果进行优化，纠正匹配错误的点，去除对结果影响极小的冗余匹配关系，输出高精度的匹配结果。

实施例一

本实施例公开了一种视网膜半自动分割方法，如图1所示，包括以下步骤：

步骤1：选取一张有m行n列像素的视网膜图像，将图像中每一列的所有像素定义为f_α＝[x₁,x₂,x₃,…,x_m]^T，α∈{1,n}，x_i∈{0,255}。

步骤2：对选取的视网膜图像进行预处理；

步骤2.1：设置5×5个像素大小的正方形滑动窗口，对图像进行均值滤波，该窗口从图像左上角开始滑动，每次移动一个像素位置，计算窗口内25个像素值的均值，用该均值替换窗口中心点像素值，该步骤循环进行直到窗口经过图像上所有点。

步骤2.2：对滤波之后的图像进行对比度增强，计算每一列像素灰度值的中值p_m，每一列中大于该列中值的像素值不变，小于中值的像素灰度值设置为0，图2和图3是本模型对视网膜OCT图像进行预处理的结果图，其中图2是健康视网膜OCT图像，图3为色素上皮层隆起的视网膜OCT图像。

步骤3：计算f_α中每个像素点的描述子；

步骤3.1：给定f_α＝[x₁,x₂,x₃,…,x_m]^T，对于任一元素x_i，取距离该点最近的l个点构成子序列p_i。将p_i分成s个可重叠的小区间，注意每个小区间的维度小于l。

步骤3.2：用最小二乘法在每个小区间内拟合条直线，该区间的特征标记为直线的斜率。最后，将子序列p_i内所有小区间的斜率整合到一起，得到一个s维的向量描述子，d_i＝f_slope(p_i)，表示点x_i的局部曲线特征。

步骤4：对OCT视网膜图像上所有曲线进行两两分组，利用动态时间规整算法求每一组曲线上点的对应关系，得到成对曲线之间的空间对应关系矩阵L_αβ∈R^mⅹm，其中1≤α，β≤n；图4是利用动态时间规划模型匹配两个灰度值曲线的一个实例。

步骤4.1：对于每一组曲线来说，首先利用欧氏距离计算两条曲线上各局部特征的相似度：

sim_(i,j)＝||d_x-t_y||²

其中，i,j为相似度矩阵S_αβ∈R^m×m中元素的位置，d_x和t_y分别表示来自两条曲线的描述子，x,y分别表示描述子在曲线上的位置。

步骤4.2：利用路径规划策略找到相似度矩阵S_αβ中相似度最小的路径：

其中，L^*(P_r)表示到第r步所需要的损失函数。另外

ω：{1：m}×{1:m}→{[0,1]^T,[1,0]^T,[1,1]^T}表示连续步骤间的转化方式，例如：

p_t+1＝p_t+ω(p_t)。为了使得该动态规划过程能有效进行，计算中做了以下约束：

(1)边界条件：所有规划路径都是从点(1,1)开始到(m,m)结束，即：从步骤P₁＝[1,1]^T开始到步骤P_k＝[m,m]^T结束。

(2)单一性：如果r₁≤r₂，则p_r1-p_r2≤0。

(3)连续性：P_r-P_r-1∈{[0,1]^T,[1,0]^T,[1,1]^T}，图5展示了动态时间规划模型进行曲线匹配时规划路径示意图。

步骤4.3：在满足步骤4.2所述约束的条件下实现目标函数：

其中，P＝[p^α,p^β]^T∈R^2×k表示两个描述子序列索引的对应关系，k表示对齐两个描述子序列所需要的步数。当D_α中的第i个描述子与D_β中的第j个描述子在第t步对应时，则

步骤4.4：通过上述过程得到的曲线上点的索引对应关系矩阵P，计算得到对应矩阵L_αβ∈R^m×m，当且仅当存在时，L_αβ中第(i,j)元素为1。把每一组曲线的对应矩阵计算出来后存储到一个连接矩阵L∈{0,1}^mn×mn中，其中，当α＝β时，L_αβ＝1。

步骤5：输入曲线两两匹配对应矩阵L，通过循环一致性的约束进行多曲线上点的匹配优化，去除L中匹配错误以及没有匹配意义的对应关系，输出精确度更高的对应矩阵X∈{0,1}^mn×mn。

步骤5.1：用公式：将X_αβ和L_αβ之间的一致性最大问题作为多个线性分配问题来解决，保证矩阵恢复无限接近真实的对应关系；利用α＜1,X＞矩阵的稀疏性；用核范数||X||_*最小化X的秩。由此构建优化目标函数：

这里，＜.,.＞表示内积，α是平衡稀疏结构成都的正则项参数，λ为核函数的权重。其中，X需满足以下约束：(1)自身匹配被记为单位矩阵。(2)X必须是对称矩阵。(3)X中的值必须在[0,1]之间。

步骤5.2：为了表示简单，我们用Z表示α1-L，因此我们优化函数简化为：

其中表示满足上述三个约束条件的所有矩阵的集合。

步骤5.3：改变变量形式令X＝AB^T，其中A,B∈R^N×u，u＜N。在这样的条件下，以下等式成立：

结合公式(4)和公式(5)目标函数变形为：

步骤5.4：用交替方向乘子算法(ADMM)来解决公式(6)最小化问题。

步骤6：手动指定点的位置通过坐标形式给出，所述指定点为边界点，根据指定点坐标位置和点的空间对应关系将与该点对应一致的点找到，用三次B样条插值进行拟合，得到平滑的分割线。

上述步骤所述的OCT眼底图像分割结果如图6和图7所示，其中图6为色素上皮层隆起病变图像分割结果，图7为上皮层下积液分割结果。

实施例二

本实施例的目的是提供一种计算装置。

一种基于曲线群匹配的OCT眼底图像半自动分割装置，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，所述处理器执行所述程序时实现以下步骤，包括：

接收一幅m行n列的眼底图像并进行预处理；

确定所述眼底图像每一列的灰度值向量曲线f_α＝[x₁,x₂,x₃,…,x_m]^T，α∈{1,n}，x_i∈{0,255}，计算灰度值向量曲线f_α中每个像素点的描述子，得到所述眼底图像的描述子矩阵；

基于所述描述子矩阵，通过动态时间规划算法将所有灰度值向量曲线进行两两匹配，得到成对曲线像素点之间的空间对应关系矩阵；

接收用户手动指定的像素点，根据所述点的坐标和所述空间对应关系矩阵，找到与所述点对应的所有像素点，拟合得到平滑的分割线。

实施例三

本实施例的目的是提供一种计算机可读存储介质。

一种计算机可读存储介质，其上存储有计算机程序，该程序被处理器执行时执行以下步骤：

接收一幅m行n列的眼底图像并进行预处理；

确定所述眼底图像每一列的灰度值向量曲线f_α＝[x₁,x₂,x₃,…,x_m]^T，α∈{1,n}，x_i∈{0,255}，计算灰度值向量曲线f_α中每个像素点的描述子，得到所述眼底图像的描述子矩阵；

基于所述描述子矩阵，通过动态时间规划算法将所有灰度值向量曲线进行两两匹配，得到成对曲线像素点之间的空间对应关系矩阵；

接收用户手动指定的像素点，根据所述点的坐标和所述空间对应关系矩阵，找到与所述点对应的所有像素点，拟合得到平滑的分割线。

以上实施例二和三的装置中涉及的各步骤与方法实施例一相对应，具体实施方式可参见实施例一的相关说明部分。术语“计算机可读存储介质”应该理解为包括一个或多个指令集的单个介质或多个介质；还应当被理解为包括任何介质，所述任何介质能够存储、编码或承载用于由处理器执行的指令集并使处理器执行本发明中的任一方法。

本发明的技术方案保证了OCT眼底图像分割的准确性，具有较强的自适应性，能够分割视网膜上明显层状结构的同时也能够处理中央凹处的重叠区域，并且可以对病变区域进行准确分割。

本领域技术人员应该明白，上述本发明的各模块或各步骤可以用通用的计算机装置来实现，可选地，它们可以用计算装置可执行的程序代码来实现，从而，可以将它们存储在存储装置中由计算装置来执行，或者将它们分别制作成各个集成电路模块，或者将它们中的多个模块或步骤制作成单个集成电路模块来实现。本发明不限制于任何特定的硬件和软件的结合。

上述虽然结合附图对本发明的具体实施方式进行了描述，但并非对本发明保护范围的限制，所属领域技术人员应该明白，在本发明的技术方案的基础上，本领域技术人员不需要付出创造性劳动即可做出的各种修改或变形仍在本发明的保护范围以内。