眼底OCT图像的视神经纤维层自动分割方法

摘要

本发明公开了一种眼底OCT图像的视神经纤维层自动分割方法，包括以下步骤：A.确定视网膜区域的内层轮廓u(s)＝[x(s)，y(s)]和外层轮廓v(s)＝[x(s)，y(s)]；B.以局部结构数据p(x，y)的方式描述出视网膜区域，并在构建水平集函数φ的基础上，建立基于局部结构数据p(x，y)的能量泛函(公式见上式)；C.最小化能量泛函，求解出φ(x，y)＝0时的分界曲线m＝{(x，y)|φ(x，y)＝0}，将分界曲线m(s)与外层轮廓v(s)之间的区域确定为视神经纤维层区域。本发明提供的方法能对眼底OCT图像的视神经纤维层进行可靠、精确的分割。

说明书

【技术领域】

本发明涉及图像处理，特别是一种眼底OCT图像的视神经纤维层自动分割方法。

【背景技术】

光学相干层析技术(Optical Coherence Tomography，OCT)是近年来迅速发展起来的一种成像技术，它在保证一定探测深度的同时，达到了高于其他方法的分辨率和灵敏度，从而以其高分辨率，高灵敏度，非侵入，非接触等特点，在诸多成像方法中占据了一席之地，并有着很好的发展前景，也受到了人们越来越多的关注。OCT利用宽带光源的弱相干特性，通过测量样品后向散射光的干涉信号，对生物组织的内部微观结构进行高分辨率层析成像。与传统的光学显微镜相比，OCT利用宽带光的弱相干性来抑制背景能量，可获得比传统的超声成像高1～2个数量级的分辨率(μm级)，以及很高的信噪比(超过100dB)。OCT影像分析方法的研究已是医学影像分析中的一个热点研究领域。OCT可用于诊断诸如青光眼、糖尿病水肿等需要定量测量视网膜变化的疾病。在青光眼的情况下，该技术可使医师能够掌握视神经纤维层的变化情形，而不必再去测量眼压及视场区域的变化。在糖尿病的情形下，眼科医生可以对视网膜的肿胀进行定量测量。眼底OCT影像分析中的视神经纤维层(retinal nerve fibre layer，RNFL)厚度计算分析在眼科疾病诊断中有着重要的作用，RNFL的厚度异常与许多眼部疾病有着直接的关系。如何精确可靠计算RNFL是OCT图像分析中的一个难点问题，主要表现在两点：1)本身RNFL的边界非常微弱而且容易与其它边界混淆，很难通过简单的方法找到RNFL；2)根据解剖学的原理RNFL应该具有双峰结构，通常针对亮度图像很难分割出具有双峰结构的RNFL。

大多数OCT系统后端分析所针对的OCT图像数据都简化为反映干涉能量的标量图像如前所述OCT图像反映的是空间光的相干能量信息，所以斑点噪声在图像重建时是无法避免的，因而降低了对生物组织描述的精度和数据处理得可靠性。一种方法从光学的角度出发，将生物组织看成大量的散射物质，最终的成像数据y(n)可以看作是散射信号x(n)与点扩散函数h(n)卷积而成，如式(1)所示。

$y\left(n\right)=x\left(n\right)\otimes h\left(n\right)---\left(1\right)$

这里点扩散函数h(n)同时可以描述为一高斯序列g(n)和一复正弦信号的乘积，如式(2)所示。

h(n)＝g(n)exp(j2πn/τ)    (2)

其中τ为光源的波长。

另一种方法将OCT图像噪声近似为乘性和加性噪声的共同作用，可以描述为式(3)。

f(x，y)＝s(x，y)n₁(x，y)+n₂(x，y)    (3)

这里s(x，y)为理想的OCT数据，f(x，y)为观测值，n₁(x，y)，n₂(x，y)分别是乘性和加性噪声，(x，y)为空间坐标。

OCT数据分析如分割、轮廓提取等，如果采用模型(2)和(3)会很复杂，实际应用时并不现实。现大多数OCT图像分析都是针对反映干涉能量的标量图像，同时利用传统的灰度图像处理技术。使用这种简化了的灰度图像在处理、分析时会降低对生物组织固有特征的描述能力，同时不能很好解决数据的噪声问题。

【发明内容】

本发明的主要目的就是解决现有技术中的问题，提供一种能对眼底OCT图像的RNFL进行可靠、精确分割的自动分割方法。

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

一种眼底OCT图像的视神经纤维层自动分割方法，包括以下步骤：

A、确定视网膜区域的内层轮廓u(s)＝[x(s)，y(s)]和外层轮廓v(s)＝[x(s)，y(s)]；

B、以局部结构数据p(x，y)的方式描述出视网膜区域，并在构建水平集函数φ的基础上，建立基于局部结构数据p(x，y)的能量泛函：

其中，d(·)为描述不同局部结构数据之间的差异的距离函数，为视神经纤维层区域Ω₁内的局部结构数据的均值，为内外层轮廓之间的非视神经纤维层区域Ω₂内的局部结构数据的均值，φ为定义在图像大小范围内的连续函数，H为指示性函数，H(x)＝1，x≥0and H(x)＝0，x＜0；

C、最小化能量泛函，求解出φ(x，y)＝0时的分界曲线将分界曲线m(s)与外层轮廓v(s)之间的区域确定为视神经纤维层区域。

优选地，所述步骤C中，通过以下扩散方程来求解能量泛函：

其中δ为Dirac函数。

优选地，所述局部结构数据为2维图像局部结构基于二阶张量来描述的张量数据：

$T=G\otimes {(▿I)}^T(▿I)=G\otimes \left(\begin{array}{cc}{I}_x^2& I_x{I}_y\\ I_x{I}_y& I_y^2\end{array}\right)$

其中I为原始OCT数据，G为高斯核函数，张量T为一个对称2×2半正定矩阵。

本发明的有益效果是：

由于OCT图像任意坐标下像素亮度难以精确反映该位置下生物组织特性，同时，OCT图像数据其局部结构特性相对于单个像素亮度对于分析、计算则更加可靠，而信息更丰富，因此，本发明使用数据局部结构信息作为OCT图像分析和处理的对象，并引入水平集(level set)的策略，将概率函数作为新的OCT数据描述手段，结合水平集策略对RNFL进行分割，能够达到从视网膜OCT图像中可靠、精确地分割和提取视神经纤维层轮廓的效果。

【附图说明】

图1为视神经纤维层在视网膜中的分布示意图；

图2为本发明一种实施例的视神经纤维层自动分割方法的流程图；

图3为本发明一种实施例中采用的张量数据的几何意义示意图；

图4为利用张量数据描述OCT图像的分析示意图；

图5a-图5b展示了使用本发明自动分割方法得到视神经层厚度的三例结果。

【具体实施方式】

本发明的特征及优点将通过实施例结合附图进行详细说明。

请参考图2，一种实施例的视神经纤维层自动分割方法的流程如下。

假设已经知道视网膜区域内外层轮廓，分别表示为：u(s)＝[x(s)，y(s)]， v(s)＝[x(s)，y(s)]，RNFL的分界曲线为m(s)＝[x(s)，y(s)]，u，v之间的区域用Ω表示，RNFL的区域用Ω₁表示，非RNFL的区域用Ω₂表示。用p(x，y)来表示OCT图像局部结构数据。

根据水平集理论建立基于非欧数据p(x，y)能量范函：

其中，d(·)为距离函数，描述不同局部结构数据之间的差异，为Ω₁内局部结构数据的均值，为Ω₂内局部结构数据的均值，φ为水平集函数，为定义在图像大小范围内的连续函数，H为指示性函数，H(x)＝1，x≥0and H(x)＝0，x＜0。

根据Euler-Lagrange方程，方程(4)的解可以通过以下扩散方程(5)实现。

δ为Dirac函数。

方程(5)解可以通过线性迭代实现，通过一定数量的迭代之后，φ(x，y)为0的曲线即为我们要寻找的分界曲线为m(s)＝[x(s)，y(s)]，即m(s)与v(s)之间的部分即为RNFL区域。

在不同是实施例中，可采用多种描述图像局部结构的方法：如基于局部方向特性、基于局部轮廓分布特性、以及基于局部图像blob拓扑性质、基于纹理特性、基于二阶张量等描述手段。其中，基于二阶张量的描述方法由于其可以采用线性数据量化描述，相对其它特征描述更便于计算和分析，因此，在优选的实例中，将OCT图像描述为张量场，这样，式(4)中的p(x，y)为局部张量数据，d(·)为描述不同局部张量数据之间差异的距离函数，可用矩阵的Frobenius范数来描述不同张量数据之间的差异，为Ω₁内局部张量数据的均值，表示Ω₂内局部张量数据的均值。

数学上，张量是由欧式空间的切矢量与协变切矢量作用产生的，对于图像结构张量的构造同样遵循上面的原则。假设2D图像点x，y)的切矢量为$▿=({\partial }_x,{\partial }_y),$按照张量的定义2D图像的局部张量为：

$T=G\otimes {(▿I)}^T(▿I)=G\otimes \left(\begin{array}{cc}{I}_x^2& I_x{I}_y\\ I_x{I}_y& I_y^2\end{array}\right)---\left(6\right)$

这里I为原始OCT数据，G为高斯核函数。张量T包含了某一空间位置的局部结构特征。

如图3和图4所示，张量T为一个对称2×2半正定矩阵，其几何含义可以描述为空间的一个椭圆，主轴长度对应T的特征值，相应的特征向量对应主轴方向。根据T的几何含义可以将T写为：$T={\lambda }_1{e}_1^T{e}_1+{\lambda }_2{e}_2^T{e}_2,$λ₁，λ₂和e₁，e₂为相应的特征值和特征矢量。

利用张量描述时，OCT图像分析可以看作是在每个坐标位置下对λ₁，λ₂和e₁，e₂的计算和处理。

图5a-图5b展示了三例使用本发明自动分割方法得到的视神经层厚度的结果。其中，横坐标表示视神经层的扫描位置，纵坐标表示相应位置上的视神经层厚度，位于上方的阴影区域表示正常的视神经层厚度范围，黑色曲线表示采用本发明的方法计算出的视神经层厚度值，图中可见，黑色曲线的纵坐标值即得到的厚度值大体精确地落入正常的视神经层厚度范围。

以上内容是结合具体的优选实施方式对本发明所作的进一步详细说明，不能认定本发明的具体实施只局限于这些说明。对于本发明所属技术领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干简单推演或替换，都应当视为属于本发明的保护范围。