CN102393956B - 基于径向对称性变换的视网膜图像滤波方法 - Google Patents

Abstract

本发明涉及一种基于径向对称性变换的视网膜图像滤波方法，包括步骤：基于海森矩阵的多尺度滤波方法得到滤波图像、血管的尺寸和血管的方向；Canny边缘检测算法得到边缘图像，计算梯度方向；再计算径向对称性计数值和径向对称性贡献值；最后，计算径向对称性变换结果，达到滤波效果。本发明方法能有效地滤除视网膜图像中的绝大部分非血管结构，有利于视网膜图像的后续处理。在视网膜图像血管分割、脑血管图像分割以及心脏图像血管分割等领域有着重要的应用价值。

Description

基于径向对称性变换的视网膜图像滤波方法

技术领域

本发明涉及图像处理技术领域，特别涉及一种基于径向对称性变换的视网膜图像滤波方法。

背景技术

在人类从外界接受的各种信息中，约有80％以上的信息来自于视觉。据统计我国每年大约有170万人因各种眼部疾病导致失眠，而这其中的相当多数可以通过早期的诊断与治疗来预防。在临床中，视网膜图像被广泛应用于诊断眼部的疾病，包括青光眼、糖尿病视网膜病变(DR)、老年黄斑变性(ARMD)等。视网膜图像通过眼底相机进行二维成像，由于现有的影像仪器由于固有的或设计上的原因，都存在着若干缺陷，表现在视网膜图像上有以下几个方面：图像的退化、图像发生几何畸变、图像发生灰度畸变等。

发明内容

(一)要解决的技术问题

视网膜图像中存在很多病理改变和其他非血管结构(如视盘)，使得现有的方法很难准确地提取出视网膜图像中的血管结构，本发明的目的在于提供一种基于径向对称性变换的视网膜图像滤波方法，使得视网膜图像中的非血管结构引起的噪声能够得到抑制，以便于血管分割方法能够从视网膜图像中准确地提取出血管。

(二)技术方案

为达到上述目的，本发明提供了一种基于径向对称性变换视网膜图像滤波方法，该方法包括以下步骤：

步骤S1：用基于海森矩阵的多尺度滤波方法对视网膜灰度图像I(x，y)进行计算滤波，得到滤波图像V(x，y)、血管半径r(x，y)和血管的方向

x，y表示坐标；

步骤S2：用Canny边缘检测算法对视网膜灰度图像I(x，y)进行检测，得到边缘图像e(x，y)，计算并得到梯度方向

步骤S3：利用梯度方向

计算并得到径向对称性计数图像O(x，y)；

步骤S4：利用滤波图像V(x，y)、血管半径r(x，y)、血管方向和径向对称性计数图像O(x，y)计算并得到径向对称性贡献值F(x，y)；

步骤S5：利用血管方向

并从径向对称性贡献值F(x，y)选取方向对应的角度A(x，y)所在区间[d，d+Δd)内的像素作为方向性图像F_d(x，y)，计算并得到径向对称性变换结果S(x，y)，其中，d表示角度，取值范围为d＝0°，Δd°，2Δd°，…，167.14°；Δd表示角度步长。

优选地，所述视网膜灰度图像I(x，y)的滤波图像V(x，y)表示为：

$V(x,y)=\underset{s_{\min }\≤s\≤s_{\max }}{\max }\left[s^{2}v(x,y,s)\right]$

其中，s表示尺度或方差；s_min，s_max分别表示尺度空间的最大和最小尺度，滤波图像V(x，y)公式具有尺度选择特性：当血管半径与尺度匹配时，滤波数值最大；v(x，y，s)表示尺度s下的滤波图像。

优选地，所述尺度s下的滤波图像v(x，y，s)表示为：

其中，λ₁(x，y，s)，λ₂(x，y，s)为尺度s下的海森矩阵两个特征值(|λ₁(x，y，s)|≥|λ₂(x，y，s)|)；

为尺度s下的海森矩阵的两个特征向量；c为一常数。

优选地，所述尺度s下的海森矩阵H(x，y，s)表示为：

$H(x,y,s)=\left[\begin{array}{cc}{I}_{\mathrm{xx}}(x,y,s)& I_{\mathrm{xy}}(x,y,s)\\ I_{\mathrm{xy}}(x,y,s)& I_{\mathrm{yy}}(x,y,s)\end{array}\right]$

其中，

I_xx(x，y，s)表示尺度s下的二次横向图像， $I_{\mathrm{xx}}(x,y,s)=\frac{\∂G(x,y,s)}{\∂x^2}*I(x,y);$

I_xy(x，y，s)表示尺度s下的二次交叉图像， $I_{\mathrm{xy}}(x,y,s)=\frac{\∂G(x,y,s)}{\∂\mathrm{xy}}*I(x,y);$

I_yy(x，y，s)表示尺度s下的二次纵向图像， $I_{\mathrm{yy}}(x,y,s)=\frac{\∂G(x,y,s)}{\∂y^2}*I(x,y);$

$G(x,y,s)=\frac{1}{2{\mathrm{\πs}}^2}{e}^{-\frac{x^2+y^2}{2s^2}};$ *表示卷积。

优选地，所述梯度方向

表示为：

$\stackrel{\→}{g}(x,y)=\stackrel{\→}{G}(x,y)/\left|\right|\stackrel{\→}{G}(x,y)\left|\right|$

其中，为视网膜灰度图像I(x，y)的梯度场；||·||表示计算幅值。

优选地，所述计算并得到径向对称性计数图像O(x，y)的步骤包括：

步骤S31：构造一幅和视网膜灰度图像I(x，y)尺寸大小一样的计数图像O(x，y)，初始化为0；

步骤S32：根据血管半径r(x，y)和梯度方向

计算边缘图像e(x，y)像素所影响到像素的坐标：

其中，

表示被影响到的像素坐标；

表示下取整；

m＝0，Δm，2Δm，…，2r(x，y)；Δm是步长；上式计算得到一次坐标，则计数图像O(x，y)中相应像素的计数值累加1。

优选地，所述径向对称性贡献值F(x，y)表示为：

$F(x,y)=V(x,y)\×{\left(\frac{O(x,y)}{k}\right)}^{\mathrm{\γ}}$

其中，k，γ是一个常数。

优选地，所述径向对称性变换结果S(x，y)表示为：

其中，L_n(d)是线状方向性高斯核；F_d(x，y)是方向性图像，其尺寸和F(x，y)一致，其像素值初始化为0，使用血管方向

计算角度图像4(x，y)，若A(x，y)在区间[d，d+Δd)内，则把F(x，y)对应的像素值传递给F_d(x，y)；n表示线状方向性高斯核的大小。

优选地，所述线状方向性高斯核L_n(d)是一个n×n大小的矩阵，该矩阵只有一条直线上的值为非零，且其值按高斯函数分布；以一个7×7矩阵且是0°方向为例，第4行元素的值按高斯分布，第4行第4列元素的值最大，矩阵中其它元素全为0，则得到14个7×7线性方向性高斯核L₇(0°)，L₇(12.86°)，L₇(25.72°)，…，L₇(167.14°)。

(三)有益效果

1、本发明提供的这种基于径向对称性变换的视网膜图像滤波方法，是视网膜血管提取的预处理技术，能大大地降低血管提取的难度，提高视网膜血管提取的准确度。

2、本发明提供的这种基于径向对称性变换的视网膜图像滤波方法，能自动地滤除视网膜图像中的非血管结构引起的噪声。在视网膜图像血管分割、脑血管图像分割和心脏图像血管分割等领域有着重要的应用价值。

本发明旨在滤除视网膜图像中的噪声，以便于准确地分割视网膜血管。本方法结合基于海森矩阵的多尺度滤波方法，同时利用径向对称性变换进一步地滤除多尺度过程带来的部分噪声，以便得到清晰的血管网络。实验结果显示，本方法可以抑制视网膜图像中病理改变和视盘等引起的噪声。

附图说明

图1为本发明提供的基于径向对称性变换的视网膜图像滤波方法的流程图。

图2a至图2c为本发明提供的径向对称性变换的示意图。

图3a至图3d为本发明提供的基于海森矩阵的多尺度滤波方法的实例。

图4a至图4d为本发明提供的基于径向对称性变换的滤波实例。

具体实施方式

下面将结合附图对本发明加以详细说明：

本发明的核心思想是提出一种基于径向对称性变换的视网膜图像滤波方法。图1为本发明提供的基于径向对称性变换的视网膜图像滤波方法的流程图，该方法基本思路为：用基于海森(Hessian)矩阵的多尺度滤波方法检测血管得到滤波图像V(x，y)，血管的尺寸r(x，y)和血管的方向

检测出视网膜图像中的边缘e(x，y)和计算梯度方向

结合血管的尺寸、边缘和梯度方向等信息，计算径向对称性计数图像O(x，y)；再利用多尺度滤波图像和径向对称性计数图像，计算径向对称性贡献值F(x，y)；最后，根据血管的方向、线状方向性高斯核和径向对称性贡献值，可计算出径向对称性变换结果S(x，y)。

以下结合具体的实施例对根据本发明提供的这种基于径向对称性变换的视网膜图像滤波方法进行详细描述。

步骤S1：用基于海森矩阵的多尺度滤波方法得到滤波图像V(x，y)、血管的尺寸r(x，y)和血管的方向

x，y表示坐标。

基本思路：首先，对视网膜灰度图像使用多尺度表示可得到图像中每个像素的海森矩阵；然后，计算海森矩阵的特征值和特征向量；最后，构建具有尺度选择性的响应函数检测出血管结构及其半径，滤波前后的图像见图3a至图3d。

一幅视网膜灰度图像I(x，y)，其多尺度表示为：

$I_{\mathrm{xx}}(x,y,s)=\frac{\∂G(x,y,s)}{\∂x^2}*I(x,y)$

$I_{\mathrm{xy}}(x,y,s)=\frac{\∂G(x,y,s)}{\∂\mathrm{xy}}*I(x,y)$

$I_{\mathrm{yy}}(x,y,s)=\frac{\∂G(x,y,s)}{\∂y^2}*I(x,y)$

其中，s表示方差或尺度；*表示卷积；I_xx(x，y，s)表示尺度s下的二次横向图像；I_xy(x，y，s)表示尺度s下的二次交叉图像；I_yy(x，y，s)表示尺度s下的二次纵向图像。可得到尺度s下的海森矩阵：

$H(x,y,s)=\left[\begin{array}{cc}{I}_{\mathrm{xx}}(x,y,s)& I_{\mathrm{xy}}(x,y,s)\\ I_{\mathrm{xy}}(x,y,s)& I_{\mathrm{yy}}(x,y,s)\end{array}\right]$

并求该矩阵的特征值λ₁(x，y，s)，λ₂(x，y，s)(|λ₁(x，y，s)|≥|λ₂(x，y，s)|)和特征向量

利用特征值即可构建多尺度滤波的响应：

$V(x,y)=\underset{s_{\min }\≤s\≤s_{\max }}{\max }\left[s^{2}v(x,y,s)\right]$

其中，s_min，s_max分别表示尺度空间的最大和最小尺度，上式具有尺度选择特性：当血管的半径与尺度s匹配时，滤波数值最大；v(x，y，s)表示尺度s下的滤波图像：

其中c为一常数。多尺度选择后得到的s即血管的半径r(x，y)，

作为血管的方向

步骤S2：Canny边缘检测算法得到边缘图像e(x，y)，计算梯度方向

用Canny边缘检测方法得到灰度图像I(x，y)的边缘图像e(x，y)，并计算灰度图像I(x，y)的梯度方向：

$\stackrel{\→}{g}(x,y)=\stackrel{\→}{G}(x,y)/\left|\right|\stackrel{\→}{G}(x,y)\left|\right|$

其中，为视网膜灰度图像I(x，y)的梯度场；||·||表示计算幅值。

步骤S3：计算径向对称性计数图像O(x，y)

基本思路：根据视网膜的边缘图像e(x，y)、梯度方向

和血管的半径r(x，y)，计算视网膜的边缘像素所影响的像素的坐标，找到该像素后其计数值则累计1.对于具有双边结构的血管来说，血管像素的径向对称性计数值来自于血管两边的累加，是单边结构(例如病理改变和视盘)的两倍。

首先，构造一幅和I(x，y)尺寸大小一样的计数图像O(x，y)，初始化为0。然后，上述步骤S1中的半径r(x，y)和上述步骤S2中的梯度方向计算边缘图像e(x，y)像素所影响到像素的坐标：

其中，

表示被影响到的像素坐标；表示下取整；

m＝0，Δm，2Δm，…，2r(x，y)；Δm是步长。上式计算得到一次坐标，则计数图像O(x，y)中相应像素的计数值累加1。

步骤S4：计算径向对称性贡献值F(x，y)

结合多尺度滤波图像V(x，y)和径向对称性计数图像O(x，y)，可计算出径向对称性贡献值：

$F(x,y)=V(x,y)\·{\left(\frac{O(x,y)}{k}\right)}^{\mathrm{\γ}}$

其中，k、γ是一个常数；

步骤S5：计算径向对称性变换结果S(x，y)

基本思路：由于上述步骤S4的F(x，y)是散乱的点，无法得到光滑的血管结构，因此，本发明构造一组线状方向性高斯核与径向对称性贡献值卷积，可得到光滑的血管。

利用步骤S1中计算得到的血管方向

从F(x，y)选取方向对应的角度A(x，y)所在区间[d，d+Δd)内的像素作为一副新的图像F_d(x，y)，计算最终的径向对称性变换：

其中，L_n(d)是线状方向性高斯核，F_d(x，y)是方向性图像，其尺寸和F(x，y)一致，其像素值初始化为0，使用血管方向

计算角度图像A(x，y)，若A(x，y)在区间[d，d+Δd)内，则把F(x，y)对应的像素值传递给F_d(x，y)；d表示角度，取值范围为d＝0°，Δd°，2Δd°，…，167.14°；Δd表示角度步长；n表示线状方向性高斯核的大小。图2a为本发明提供的径向对称性变换的示意图。线状方向性高斯核L_n(d)是一个n×n大小的矩阵，该矩阵只有一条直线上的值为非零，且其值按高斯函数分布。以一个7×7矩阵且是0°方向(如图2b所示)，第4行元素的值按高斯分布(第4行第4列元素的值最大)，矩阵中其它元素全为0，则可得到14个7×7线性方向性高斯核L₇(0°)，L₇(12.86°)，L₇(25.72°)，…，L₇(167.14°)(如图2c所示)。

运行结果：

为验证本发明的有效性和实用性，我们在两个国际公认的视网膜图像数据库(STARE和DRIVE)上进行了实验。这两个数据库都分别提供20幅视网膜图像用于算法测试，两个手动分割的血管网络数据集作为参考。

图3a至图3d和图4a至图4d给出了一组具体的例子。图3a为带有视盘的视网膜图像，图3b为带有病理改变的视网膜的图像，图3c为图3a使用基于海森矩阵的多尺度滤波方法得到的图像，图3d为图3b为使用基于海森矩阵的多尺度滤波方法得到的图像。图4a为基于径向对称性变换的滤波图像，图4b为基于径向对称性变换的滤波图像，图4c为手动分割结果，图4d为手动分割结果。在本实验中，本滤波方法都能有效地把绝大多数的噪声滤除，特别地，包括病理改变和视盘引起的噪声。从原图和滤波后图像的对比可以看出，经过本发明的处理，血管网络已经突显出来。

实验表明，本发明方法-径向对称性变换的视网膜图像滤波方法-抑制病理改变和视盘引起的噪声，达到了滤除视网膜中绝大部分噪声的目的。

以上所述，仅为本发明中的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉该方法的人在本发明所揭露的方法范围内，可理解想到的变换或替换，都应涵盖在本发明的包含范围之内。

Claims (8)

1.一种基于径向对称性变换的视网膜图像滤波方法，其特征在于，包括以下步骤：

步骤S1：用基于海森矩阵的多尺度滤波方法对视网膜灰度图像I(x,y)进行计算滤波，得到滤波图像V(x,y)、血管半径r(x,y)和血管的方向

x,y表示坐标；

步骤S2：用Canny边缘检测算法对视网膜灰度图像I(x,y)进行检测，得到边缘图像e(x,y)，计算并得到梯度方向

；

步骤S3：利用梯度方向

计算并得到径向对称性计数图像O(x,y)；

步骤S4：利用滤波图像V(x,y)、血管半径r(x,y)、血管方向和径向对称性计数图像O(x,y)计算并得到径向对称性贡献值F(x,y)；

步骤S5：利用血管的方向

，并从径向对称性贡献值F(x,y)选取方向对应的角度A(x,y)所在区间[d,d+Δd)内的像素作为方向性图像Fd(x,y)，计算并得到径向对称性变换结果S(x,y)，其中，d表示角度，取值范围为d=0°,Δd°,2Δd°,…,167.14°；Δd表示角度步长；

其中，所述视网膜灰度图像I(x,y)的滤波图像V(x,y)表示为：

$(x,y)=\underset{s_{\min }\≤s\≤s_{\max }}{\max }\left[s^{2}v(x,y,s)\right]$

其中，s表示尺度或方差；s_min,s_max分别表示尺度空间的最大和最小尺度，滤波图像V(x,y)公式具有尺度选择特性：当血管半径与尺度匹配时，滤波数值最大；v(x,y,s)表示尺度s下的滤波图像。

2.根据权利要求1所述的基于径向对称性变换的视网膜图像滤波方法，其特征在于，所述尺度s下的滤波图像v(x,y,s)表示为：

其中，λ₁(x,y,s),λ₂(x,y,s)为尺度s下的海森矩阵两个特征值（|λ₁(x,y,s)|≥|λ₂(x,y,s)|）；

为尺度s下的海森矩阵的两个特征向量；c为一常数。

3.根据权利要求1所述的基于径向对称性变换的视网膜图像滤波方法，其特征在于，所述尺度s下的海森矩阵H(x,y,s)表示为：

$H(x,y,s)=\left[\begin{array}{cc}{I}_{\mathrm{xx}}(x,y,s)& I_{\mathrm{xy}}(x,y,s)\\ I_{\mathrm{xy}}(x,y,s)& I_{\mathrm{yy}}(x,y,s)\end{array}\right]$

其中，

I_xx(x,y,s)表示尺度s下的二次横向图像，

I_xy(x,y,s)表示尺度s下的二次交叉图像，I_yy(x,y,s)表示尺度s下的二次纵向图像，

$G(x,y,s)=\frac{1}{2\mathrm{\π}{s}^2}{e}^{-\frac{x^2+y^2}{2s^2}};$ *表示卷积。

4.根据权利要求1所述的基于径向对称性变换的视网膜图像滤波方法，其特征在于，所述梯度方向

表示为：

$\stackrel{\→}{g}(x,y)=\stackrel{\→}{G}(x,y)/\left|\right|\stackrel{\→}{G}(x,y)\left|\right|$

其中，为视网膜灰度图像I(x,y)的梯度场；||·||表示计算幅值。

5.根据权利要求1所述的基于径向对称性变换的视网膜图像滤波方法，其特征在于，所述计算并得到径向对称性计数图像O(x,y)的步骤包括：

步骤S31:构造一幅和视网膜灰度图像I(x,y)尺寸大小一样的计数图像O(x,y)，初始化为0；

步骤S32:根据血管半径r(x,y)和梯度方向

计算边缘图像e(x,y)像素所影响到像素的坐标：

其中，

表示被影响到的像素坐标；表示下取整；

m=0,Δm,2Δm,…,2r(x,y)；Δm是步长；上式计算得到一次坐标，则计数图像O(x,y)中相应像素的计数值累加1。

6.根据权利要求1所述的基于径向对称性变换的视网膜图像滤波方法，其特征在于，所述径向对称性贡献值F(x,y)表示为：

$F(x,y)=V(x,y)\×{\left(\frac{O(x,y)}{k}\right)}^{\mathrm{\γ}}$

其中，k，γ是一个常数。

7.根据权利要求1所述的基于径向对称性变换的视网膜图像滤波方法，其特征在于，所述径向对称性变换结果S(x,y)表示为：

其中，L_n(d)是线状方向性高斯核；F_d(x,y)是方向性图像，其尺寸和F(x,y)一致，其像素值初始化为0，使用血管方向

计算角度图像A(x,y)，若A(x,y)在区间[d,d+Δd)内，则把F(x,y)对应的像素值传递给F_d(x,y)；n表示线状方向性高斯核的大小。

8.根据权利要求7所述的基于径向对称性变换的视网膜图像滤波方法，其特征在于，所述线状方向性高斯核L_n(d)是一个n×n大小的矩阵，该矩阵只有一条直线上的值为非零，且其值按高斯函数分布。

Images