CN103198321A - 视网膜图像处理方法及装置 - Google Patents

Abstract

本发明公开了一种视网膜图像处理方法，包括：计算视网膜图像的血管结构显著性函数，对视网膜图像进行增强；对增强后的视网膜图像实施非局部均值滤波，抑制增强后的视网膜图像的噪声；采用基于图像梯度对称性的方法去除所述视网膜图像中的非血管结构的像素；采用基于图割的方法提取视网膜图像中的血管结构。本发明还公开了一种视网膜图像处理装置。本发明能有效地增强视网膜图像中的血管结构信息，最终能准确地提取视网膜图像的血管结构。

Description

视网膜图像处理方法及装置

技术领域

本发明涉及图像识别技术领域，尤其涉及一种视网膜图像处理方法及装置。

背景技术

图像增强的方法是通过一定手段对原图像附加一些信息或变换数据，有选择地突出图像中感兴趣的特征或者抑制图像中某些不需要的特征，使图像与视觉响应特性相匹配。常用的图像增强方法包括平滑和锐化。平滑一般用于消除图像噪声，常用的方法有均值滤波、中值滤波等。这类方法容易引起边缘的模糊，不利于视网膜图像中的微细结构提取。锐化的目的在于突出物体的边缘轮廓，便于目标识别，常用算法有梯度算子、拉普拉斯算子。这类基于梯度的方法对噪声比较敏感，增强后的图像往往会引入奇异的噪点或者斑块，对于视网膜图像的增强效果不是很理想。

图像分割就是把图像分成若干个特定的、具有独特性质的区域并提出感兴趣目标的技术和过程。常用的图像分割方法主要有基于阈值的分割方法、基于区域的分割方法、基于边缘的分割方法以及基于特定理论的分割方法等。

基于阈值的分割方法运算效率较高、速度快，但是在许多情况下，物体和背景的对比度在图像中的各处不是一样的，这时很难用一个统一的阈值将物体与背景分开。

基于区域的分割方法的基本思想是将具有相似性质的像素集合起来构成区域。具体先对每个需要分割的区域找一个种子像素作为生长的起点，然后将种子像素周围邻域中与种子像素有相同或相似性质的像素(根据某种事先确定的生长或相似准则来判定)合并到种子像素所在的区域中。将这些新像素当作新的种子像素继续进行上面的过程，直到再没有满足条件的像素可被包括进来。这样一个区域就长成了。基于区域的分割方法的优点是计算简单，对于较均匀的连通目标有较好的分割效果。它的缺点是需要人为确定种子点，对噪声敏感，可能导致区域内有空洞。

基于边缘的分割方法的基本思想是检测灰度级或者结构具有突变的地方，对于突变导致的灰度值不连续性可以通过求导数来检测。常用的一阶微分算子有Roberts算子、Prewitt算子和Sobel算子，二阶微分算子有Laplace算子和Kirsh算子等。基于边缘的分割方法主要缺点是这些微分算子对噪声非常敏感，只适合于噪声较小且不太复杂的图像。

视网膜图像中存在复杂的血管结构以及视盘、中央凹、黄斑等非血管结构，现有的图像处理方法很难有效地增强血管结构信息，并准确地提取出血管结构。

发明内容

本发明实施例提出一种视网膜图像处理方法及装置，使得视网膜图像中的血管结构信息在视觉效果上能显著地增强，并且能有效抑制噪声，最终能准确地提取视网膜图像的血管结构。

本发明实施例提供一种视网膜图像处理方法，包括：

获取待处理的视网膜图像的各个像素的灰度值数据；

使用所有像素的灰度值数据，计算视网膜图像的血管结构显著性函数，对视网膜图像进行增强；

对增强后的视网膜图像实施非局部均值滤波，抑制增强后的视网膜图像的噪声；

对滤波后的视网膜图像，采用基于图像梯度对称性的方法去除所述视网膜图像中的非血管结构的像素；

采用基于图割的方法，提取视网膜图像的血管结构。

相应地，本发明实施例还提供一种视网膜图像处理装置，包括：

数据获取单元，用于获取待处理的视网膜图像的各个像素的灰度值数据；

图像增强单元，用于使用所有像素的灰度值数据，计算视网膜图像的血管结构显著性函数，对视网膜图像进行增强；

滤波单元，用于对增强后的视网膜图像实施非局部均值滤波，抑制增强后的视网膜图像的噪声；

梯度对称性检测单元，用于对滤波后的视网膜图像，采用基于图像梯度对称性的方法去除所述视网膜图像中的非血管结构的像素；和，

图割提取单元，用于采用基于图割的方法，提取视网膜图像的血管结构。

本发明实施例提供的视网膜图像处理方法及装置，首先计算视网膜图像的血管结构显著性函数，使得视网膜图像中的血管结构信息在视觉效果上能显著地增强；然后对增强后的视网膜图像实施非局部均值滤波，以有效地抑制增强后的视网膜图像的噪声；再采用基于图像梯度对称性的方法去除非管状结构，最后采用基于图割的方法提取视网膜图像的血管结构。本发明能有效地增强视网膜图像中的血管结构信息，最终能准确地提取视网膜图像的血管结构。

附图说明

图1是本发明提供的视网膜图像处理方法的一个实施例的流程示意图；

图2是本发明提供的待处理的视网膜图像的一个实施例的示意图；

图3是本发明提供的增强后的视网膜图像的示意图；

图4是本发明提供的滤波后的视网膜图像的示意图；

图5是本发明提供的梯度对称性的示意图；

图6是本发明提供的提取的视网膜图像的血管结构的示意图；

图7是本发明提供的视网膜图像处理装置的一个实施例的结构示意图。

具体实施方式

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

参见图1，是本发明提供的视网膜图像处理方法的一个实施例的流程示意图。所述视网膜图像处理方法，包括以下步骤S1～S5：

S1、获取待处理的视网膜图像的各个像素的灰度值数据。

待处理的视网膜图像可以是普通的数码电子照片，例如jpg、tiff等格式的照片，还可以是其他格式的图片。如图2所示，是本发明提供的待处理的视网膜图像的一个实施例的示意图。

具体实施时，计算机加载待处理的视网膜图像，并将图像转化为二维数组，存储在计算机中。二维数组的数据为像素点的灰度值数据，数组的大小等同于图像像素尺寸大小，例如，待处理的视网膜图像的原图像大小为1024列*768行，则转换后的二维数组的大小也为1024列*768行。

S2、使用所有像素的灰度值数据，计算视网膜图像的血管结构显著性函数，对视网膜图像进行增强。

所述血管结构显著性函数是具有如下数学表达形式的函数：

$\mathrm{VL}\left(X\right)=\underset{s_{\min }<s<s_{\max }}{\max }V\left(s\right)$

其中，VL(X)是血管结构显著性函数，s是血管半径，s_min是血管半径的最小值，s_max是血管半径的最大值，V(s)是尺度s下的管状结构特征函数。

所述尺度s下的管状结构特征函数是具有如下数学表达形式的函数：

$V\left(s\right)=\frac{\left|{\mathrm{\&lambda;}}_{1,s}\right|}{2}\&CenterDot;e^{|\mathrm{\&gamma;}-\frac{\left|{\mathrm{\&lambda;}}_{1,s}\right|}{\sqrt{{{\mathrm{\&lambda;}}_{1,s}}^2+{{\mathrm{\&lambda;}}_{2,s}}^2}}|}+\frac{\left|{\mathrm{\&lambda;}}_{2,s}\right|}{2}\&CenterDot;e^{|\frac{\left|{\mathrm{\&lambda;}}_{2,s}\right|}{\sqrt{{{\mathrm{\&lambda;}}_{1,s}}^2+{{\mathrm{\&lambda;}}_{2,s}}^2}}-\mathrm{\&gamma;}|}$

其中，γ是一个经验常数(根据实际经验确定的常数)，λ_1，s和λ_2，s是视网膜图像在尺度s下的汉森(Hessian)矩阵H(X，s)的特征值，并且满足以下条件：|λ_1，s|≤|λ_2，s|。

需要说明的是，汉森(Hessian)矩阵为本技术领域的现有技术，在此不进行详细描述。关于汉森矩阵的详细资料，可参考以下文献：A.F.Frangi，W.J.Niessen，K.L.Vincken et al，“Multiscale vessel enhancement filtering”，MedicalImaging Computing and Computer Assisted Interventions，Lecture Notes inComputer Sciences，1496，pp.130-7，(1998).

在步骤S2中，计算视网膜图像的血管结构显著性函数时，是将步骤S1的二维数组中的每一个像素的灰度值数据，逐个地分别代入血管结构显著性函数中进行计算，最后计算出的“视网膜图像的血管结构显著性函数值”是一个矩阵，包含所有像素点的显著性函数值。经计算得出视网膜图像的血管结构显著性函数值，使得视网膜图像中的血管结构信息在视觉效果上能显著地增强，如图3所示。

S3、对增强后的视网膜图像实施非局部均值滤波，抑制增强后的视网膜图像的噪声。

所述非局部均值滤波是按如下数学表达式实施：

$\mathrm{NL}\left(X_i\right)=\underset{j\&Element;I}{\mathrm{\&Sigma;}}w(i,j)\&CenterDot;\mathrm{VL}\left(X_j\right)$

其中，i和j是图像下标(i，j分别代表两个不同像素的下标)，I是滤波的邻域范围，NL(X_i)是像素X_i非局部均值滤波增强后的像素值。w(i，j)是滤波权重因子。

所述滤波权重因子w(i，j)是按如下数学表达式计算：

$w(i,j)=\frac{1}{Z\left(i\right)}\&CenterDot;\exp (-\frac{{\left|\right|v\left(N_i\right)-v\left(N_j\right)\left|\right|}^2}{h^2})$

$Z\left(i\right)=\underset{j}{\mathrm{\&Sigma;}}\exp (-\frac{{\left|\right|v\left(N_i\right)-v\left(N_j\right)\left|\right|}^2}{h^2})$

其中，v(N_i)是以像素X_i为中心的邻域内像素值构成的多维特征向量，h是一个滤波经验常数。最后得到滤波后的血管结构显著性函数值，如图4所示。

在一个实施方式中，I是5×5滤波的邻域范围，v(N_i)是以像素X_i为中心的3×3邻域内像素值构成的多维特征向量。

5×5的邻域，是指以当前处理的像素(矩阵坐标为[m，n])为中心，最近邻的25个像素点(包含当前处理的像素)所构成的一个5×5子矩阵，该5×5子矩阵中的像素点的坐标[p，q]满足以下条件：m-2≤p≤m+2，n-2≤q≤n+2。3×3的邻域也类似。

S4、对滤波后的视网膜图像，采用基于图像梯度对称性的方法去除所述视网膜图像中的非血管结构的像素。

所述基于图像梯度对称性的方法去除非管状结构，包括：计算所述视网膜图像的血管对称性显著函数VL_Sym(X)；当VL_Sym(X_i)＝0时，则像素X_i不属于血管结构，去除所述像素X_i；VL_Sym(X)是按如下数学表达式计算：

VL_Sym(X)＝NL(X)·Flag_G(X)

其中，Flag_G(X)是一个梯度对称性示信函数；若像素X_i具有梯度对称性，则Flag_G(X_i)的值为1，否则Flag_G(X_i)的值为0。

X表示的是一个矩阵，是滤波后的视网膜图像中的所有像素点；而像素X_i是矩阵X中的一个像素。

所述像素X_i具有梯度对称性是指对于像素X_i，搜索其沿梯度方向的邻域[X_i-Δp，X_i+Δp]，存在满足如下数学表达式的像素X₁和X₂：

其中，Δp＝2·s(X_i)·g(X_i)，g(X_i)是像素X_i的归一化梯度向量。如图5所示。

步骤S4是对视网膜图像的血管结构的粗略分割，还不能将非血管结构的像素完全去除，仍然有部分非血管结构的像素被误认为是血管像素。需要在步骤S5中对视网膜图像的血管结构作进一步的精细分割，以准确地提取视网膜图像的血管结构。

S5、采用基于图割的方法，提取视网膜图像的血管结构。

所述采用基于图割的方法，提取视网膜图像中的血管结构，包括：

将视网膜图像分成血管像素T和非血管像素S，采用最大流/最小割的方法使得以下表达式的能量函数达到最小值：

E(L)＝E_d(L)+kE_s(L)

其中，L＝{S，T}是一个对视网膜图像的划分集合，

是将像素X_i划分为S或者T的能量代价，k是一个经验常数，

是将两个相邻像素p、q划分为不同类像素的能量代价。

所述将像素X_i划分为S或者T的能量代价，是按如下数学表达式计算：

其中，T_H是高阈值，T_L是低阈值，m是一个经验常数，μ(T)是所有血管像素的均值，μ(S)是所有非血管像素的均值。

需要说明的是，“最大流/最小割的方法”为本技术领域的现有技术，在此不进行详细描述。关于“最大流/最小割的方法”的详细资料，可参考以下文献：Y.Boykov and V.Kolmogorov，“An experimental comparison of mincut/max-flowalgorithms for energy minimization in vision”，Pattern Analysis and MachineIntelligence，IEEE Transactions on，vol.26，no.9，pp.1124-1137，2004.

所述将两个相邻像素p、q划分为不同类像素的能量代价，是按如下数学表达式计算：

$V_{\mathrm{pq}}(L_p\&NotEqual;L_q)=\frac{1}{|V\left(p\right)-V\left(q\right)|+\mathrm{\&eta;}}$

其中，V(p)，V(q)是像素值，η是一个正实数。

η是为了防止除0而添加的一个非常小的正实数，η＜＝10^-6。例如，η取值范围为[0.0000001，0.000001]。进行图像分割后，所提取的视网膜图像的血管结构如图6所示。

本发明还提供一种视网膜图像处理装置，能够实现上述实施例中的视网膜图像处理方法的所有流程。

参见图7，是本发明提供的视网膜图像处理装置的一个实施例的结构示意图。所述视网膜图像处理装置包括：

数据获取单元101，用于获取待处理的视网膜图像的各个像素的灰度值数据；

图像增强单元102，用于使用所有像素的灰度值数据，计算视网膜图像的血管结构显著性函数，对视网膜图像进行增强；

滤波单元103，用于对增强后的视网膜图像实施非局部均值滤波，抑制增强后的视网膜图像的噪声；

梯度对称性检测单元104，用于对滤波后的视网膜图像，采用基于图像梯度对称性的方法去除所述视网膜图像中的非血管结构的像素；和，

图割提取单元105，用于采用基于图割的方法，提取视网膜图像的血管结构。

本发明实施例提供的视网膜图像处理方法及装置，首先计算视网膜图像的血管结构显著性函数，使得视网膜图像中的血管结构信息在视觉效果上能显著地增强；然后对增强后的视网膜图像实施非局部均值滤波，以有效地抑制增强后的视网膜图像的噪声；再采用基于图像梯度对称性的方法去除非管状结构，最后采用基于图割的方法提取视网膜图像的血管结构。本发明能有效地增强视网膜图像中的血管结构信息，最终能准确地提取视网膜图像的血管结构。

本领域普通技术人员可以理解实现上述实施例方法中的全部或部分流程，是可以通过计算机程序来指令相关的硬件来完成，所述的程序可存储于一计算机可读取存储介质中，该程序在执行时，可包括如上述各方法的实施例的流程。其中，所述的存储介质可为磁碟、光盘、只读存储记忆体(Read-Only Memory，ROM)或随机存储记忆体(Random Access Memory，RAM)等。

以上所述是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也视为本发明的保护范围。

Claims (10)

1.一种视网膜图像处理方法，其特征在于，包括：

获取待处理的视网膜图像的各个像素的灰度值数据；

使用所有像素的灰度值数据，计算视网膜图像的血管结构显著性函数，对视网膜图像进行增强；

对增强后的视网膜图像实施非局部均值滤波，抑制增强后的视网膜图像的噪声；

对滤波后的视网膜图像，采用基于图像梯度对称性的方法去除所述视网膜图像中的非血管结构的像素；

采用基于图割的方法，提取视网膜图像的血管结构。

2.如权利要求1所述的视网膜图像处理方法，其特征在于，所述血管结构显著性函数是具有如下数学表达形式的函数：

$\mathrm{VL}\left(X\right)=\underset{s_{\min }<s<s_{\max }}{\max }V\left(s\right)$

其中，VL(X)是血管结构显著性函数，s是血管半径，s_min是血管半径的最小值，s_max是血管半径的最大值，V(s)是尺度s下的管状结构特征函数。

3.如权利要求2所述的视网膜图像处理方法，其特征在于，所述尺度s下的管状结构特征函数是具有如下数学表达形式的函数：

$V\left(s\right)=\frac{\left|{\mathrm{\&lambda;}}_{1,s}\right|}{2}\&CenterDot;e^{|\mathrm{\&gamma;}-\frac{\left|{\mathrm{\&lambda;}}_{1,s}\right|}{\sqrt{{{\mathrm{\&lambda;}}_{1,s}}^2+{{\mathrm{\&lambda;}}_{2,s}}^2}}|}+\frac{\left|{\mathrm{\&lambda;}}_{2,s}\right|}{2}\&CenterDot;e^{|\frac{\left|{\mathrm{\&lambda;}}_{2,s}\right|}{\sqrt{{{\mathrm{\&lambda;}}_{1,s}}^2+{{\mathrm{\&lambda;}}_{2,s}}^2}}-\mathrm{\&gamma;}|}$

其中，γ是一个经验常数，λ_1，s和λ_2，s是视网膜图像在尺度s下的汉森矩阵的特征值，并且满足以下条件：|λ_1，s|≤|λ_2，s|。

4.如权利要求1所述的视网膜图像处理方法，其特征在于，所述对增强后的视网膜图像实施非局部均值滤波，是按如下数学表达式实施：

$\mathrm{NL}\left(X_i\right)=\underset{j\&Element;I}{\mathrm{\&Sigma;}}w(i,j)\&CenterDot;\mathrm{VL}\left(X_j\right)$

其中，i和j是图像下标，I是滤波的邻域范围，NL(X_i)是像素X_i非局部均值滤波增强后的像素值，w(i，j)是滤波权重因子。

5.如权利要求4所述的视网膜图像处理方法，其特征在于，所述滤波权重因子w(i，j)是按如下数学表达式计算：

$w(i,j)=\frac{1}{Z\left(i\right)}\&CenterDot;\exp (-\frac{{\left|\right|v\left(N_i\right)-v\left(N_j\right)\left|\right|}^2}{h^2})$

$Z\left(i\right)=\underset{j}{\mathrm{\&Sigma;}}\exp (-\frac{{\left|\right|v\left(N_i\right)-v\left(N_j\right)\left|\right|}^2}{h^2})$

其中，v(N_i)是以像素X_i为中心的邻域内像素值构成的多维特征向量，h是一个滤波经验常数。

6.如权利要求1所述的视网膜图像处理方法，其特征在于，所述对滤波后的视网膜图像，采用基于图像梯度对称性的方法去除所述视网膜图像中的非血管结构的像素，包括：

计算所述视网膜图像的血管对称性显著函数VL_Sym(X)；当VL_Sym(X_i)＝0时，则像素X_i不属于血管结构，去除所述像素X_i；VL_Sym(X)是按如下数学表达式计算：

VL_Sym(X)＝NL(X)·Flag_G(X)

其中，Flag_G(X)是一个梯度对称性示信函数；若像素X_i具有梯度对称性，则Flag_G(X_i)的值为1，否则FlagG(X_i)的值为0；

所述像素X_i具有梯度对称性是指对于像素X_i，搜索其沿梯度方向的邻域[X_i-Δp，X_i+Δp]，存在满足如下数学表达式的像素X₁和X₂：

其中，Δp＝2·s(X_i)·g(X_i)，g(X_i)是像素X_i的归一化梯度向量。

7.如权利要求1所述的视网膜图像处理方法，其特征在于，所述采用基于图割的方法，提取视网膜图像中的血管结构，包括：

将视网膜图像分成血管像素T和非血管像素S，采用最大流/最小割的方法使得以下表达式的能量函数达到最小值：

E(L)＝E_d(L)+kE_s(L)

其中，L＝{S，T}是一个对视网膜图像的划分集合，

是将像素X_i划分为S或者T的能量代价，k是一个经验常数，

是将两个相邻像素p、q划分为不同类像素的能量代价。

8.如权利要求7所述的视网膜图像处理方法，其特征在于，所述将像素X_i划分为S或者T的能量代价，是按如下数学表达式计算：

其中，T_H是高阈值，T_L是低阈值，m是一个经验常数，μ(T)是所有血管像素的均值，μ(S)是所有非血管像素的均值。

9.如权利要求7所述的视网膜图像处理方法，其特征在于，所述将两个相邻像素p、q划分为不同类像素的能量代价，是按如下数学表达式计算：

$V_{\mathrm{pq}}(L_p\&NotEqual;L_q)=\frac{1}{|V\left(p\right)-V\left(q\right)|+\mathrm{\&eta;}}$

其中，V(p)，V(q)是像素值，η是一个正实数。

10.一种视网膜图像处理装置，包括：

数据获取单元，用于获取待处理的视网膜图像的各个像素的灰度值数据；

图像增强单元，用于使用所有像素的灰度值数据，计算视网膜图像的血管结构显著性函数，对视网膜图像进行增强；

滤波单元，用于对增强后的视网膜图像实施非局部均值滤波，抑制增强后的视网膜图像的噪声；

梯度对称性检测单元，用于对滤波后的视网膜图像，采用基于图像梯度对称性的方法去除所述视网膜图像中的非血管结构的像素；和，

图割提取单元，用于采用基于图割的方法，提取视网膜图像的血管结构。

Images