发明内容
本发明为解决上述技术问题而提供一种基于区域增长PCNN的视网膜血管分割方法。
为了实现上述发明目的,本发明采取如下技术方案:
一种基于区域增长PCNN的视网膜血管分割方法,包括:
从目标视网膜血管图像的未标记像素中选择种子点;
增大PCNN模型的连接强度,使用连接强度增大后的所述PCNN模型以所述种子点为起始点提取所述目标视网膜血管图像中的血管特征,直到增大后的连接强度大于第一预设阈值;
若本次迭代提取的血管特征不同时满足第一预设条件和第二预设条件,则将本次迭代提取的血管特征对应的像素标记为同一标签,直到所述目标视网膜血管图像的所有像素都标记有标签;
其中,所述第一预设条件为本次迭代提取的血管特征中血管边缘像素的数量占血管像素总数量的比例小于或等于第二预设阈值;
所述第二预设条件为本次迭代提取的血管特征中血管面积占所述目标视网膜血管图像面积的比例小于或等于第三预设阈值。
若本次迭代提取的血管特征不同时满足第一预设条件和第二预设条件,则将本次迭代提取的血管特征对应的像素标记为同一标签,直到所述目标视网膜血管图像的所有像素都标记有标签的步骤还包括:
若本次迭代提取的血管特征同时满足第一预设条件和第二预设条件,则将本次迭代提取的血管特征删除。
从目标视网膜血管图像的未标记像素中选择种子点的步骤之前还包括:
将待处理视网膜血管图像的各通道乘以相应预设权重后相加,获取所述待处理视网膜血管图像的灰度图像;
根据所述灰度图像,获取所述目标视网膜血管图像。
根据所述灰度图像,获取所述目标视网膜血管图像的步骤具体包括:
分别基于二维高斯滤波方法和二维Gabor滤波方法对所述灰度图像进行滤波;
将所述二维高斯滤波方法的滤波结果与二维Gabor滤波方法的滤波结果进行融合,获取融合图像;
将所述灰度图像与所述融合图像进行相减,将相减结果进行取反操作后,乘以预设掩膜图像,获取所述目标视网膜血管图像。
分别基于二维高斯滤波方法和二维Gabor滤波方法对所述灰度图像进行滤波的步骤之前还包括:
若基于四邻域法获知所述灰度图像中各像素为边界,则将所述灰度图像中为边界的各像素进行边缘膨胀;
基于对比度受限制的直方图均衡化方法对边缘膨胀后的所述灰度图像进行对比度增强,获取预处理后的所述灰度图像;
相应地,分别基于二维高斯滤波方法和二维Gabor滤波方法对所述灰度图像进行滤波的步骤具体包括:
分别基于二维高斯滤波方法和二维Gabor滤波方法对预处理后的所述灰度图像进行滤波。
若基于四邻域法获知所述灰度图像中各像素为边界,则将所述灰度图像中为边界的各像素进行边缘膨胀的步骤具体包括:
若所述灰度图像中各像素的四邻域中至少有一个值为0且不全为0,则获知所述灰度图像中各像素为边界;
将所述灰度图像中为边界的各像素进行边缘膨胀。
基于对比度受限制的直方图均衡化方法对边缘膨胀后的所述灰度图像进行对比度增强,获取预处理后的所述灰度图像的步骤具体包括:
将边缘膨胀后的所述灰度图像划分为多个块;
计算出各所述块的灰度直方图,并对各所述块的对灰度直方图进行裁剪;
对各所述块裁剪后的灰度直方图进行均衡化;
基于线性插值法将均衡化后灰度直方图对应的块进行连接,获取预处理后的所述灰度图像。
所述二维高斯滤波方法的第i个高斯卷积核K′i(x,y)的公式为:
K′i(x,y)=Ki(x,y)-mi;
Ki(x,y)=-exp(-u2/(2σ2));
其中,K
i(x,y)为第i个高斯核矩阵中的系数,m
i为第i个高斯核的均值,
为第i个高斯核的离散点,u、v分别为第i个高斯核在x轴和y轴上的值,σ为高斯核距离x轴坐标中心的偏移量,N={[u,v]||u|≤3σ,|v|≤L/2},L为在y轴方向上被截断的血管长度,A为N中的点数。
将所述二维高斯滤波方法的滤波结果与二维Gabor滤波方法的滤波结果进行融合,获取融合图像的步骤具体包括:
将所述二维高斯滤波方法的滤波结果和二维Gabor滤波方法的滤波结果分别乘以相应预设权重后相加,获取融合图像。
将相减结果进行取反操作后,乘以预设掩膜图像,获取所述目标视网膜血管图像的步骤具体包括:
对预设掩模图像相继进行开运算、闭运算和腐蚀操作;
将处理后的所述预设掩模图像与取反结果进行相乘,获取所述目标视网膜血管图像。
本发明一种视网膜血管分割方法,有益效果如下:
(1)从DRIVE数据库中随机抽出10幅图,通过准确度(Acc),灵敏度(Sen),特异性(Spe)的计算,反映出预处理中使用高斯滤波与Gabor滤波相结合的办法,比仅使用Gabor滤波或高斯滤波处理结果好。
(2)本发明在对视网膜血管预处理的基础上首先选择出一定的血管区域作为初始种子区域;接着,将带有快速连接机制的PCNN与种子区域增长思想相结合,通过设置PCNN中的连接强度系数最大值和停止条件,如果生长无效则会删除掉那部分伪血管,解决了以前方法中断点不能继续生长的问题,实现眼底图像中血管区域的自动生长,从而有效地提取出视网膜图像中的血管。
具体实施方式
一种基于区域增长PCNN的视网膜血管分割方法,括以下步骤,流程图如图1所示:
步骤1:从标准图像库DRIVE中选择一幅彩色眼底视网膜图像作为待处理视网膜血管图像,如图2所示;
步骤2:按比例提取红绿蓝通道图像Y=0.299R+0.587G+0.114B,如图3所示;
步骤3:用四邻域法判断是否为边界,若判断为边界,则将边缘膨胀,如图4所示;
步骤4:采用对比度受限制的直方图均衡化(CLAHE)的方法增强图像的对比度,如图5所示;
步骤5:采用二维高斯滤波的方法增强视网膜血管与背景之间的对比度,如图6所示
步骤6:采用二维Gabor滤波的方法进一步增强视网膜血管与背景之间的对比度,如图7所示;
步骤7:将二维Gabor与二维高斯匹配滤波结果相融合,如图8所示
步骤8:将按红绿蓝成相应比例提取的图像与经过匹配滤波处理融合的图像相减,如图9所示
步骤9:进行取反操作,乘以掩膜图像,掩膜图像如图10所示,即可得到预处理结果,即目标视网膜血管图像,如图11所示;
步骤10:将带有快速连接机制的PCNN与种子区域增长思想相结合,通过阈值操作首先选择出可靠的血管区域作为种子,然后通过PCNN中的连接强度系数和停止条件,实现血管区域的自动生长,从而有效提取视网膜血管,如图12所示。
进一步的,所述步骤2具体为:考虑到虽然绿色通道内含有大量有用信息,但是红色通道和蓝色通道仍然含有有用信息,因此从三个不同通道内按照红色提取29.9%,绿色提取58.7%,蓝色提取11.4%
进一步的,所述步骤3具体为:
步骤3.1:用四邻域的方法判断,如若四个值里有一个值为0但不全为0即为图像的边界
步骤3.2:若为边界则使用膨胀腐蚀的方法使边界处理更加平滑
进一步的,所述步骤4具体为:
步骤4.1:将原始图像分为M×N个大小相等的连续但不重叠的块;
步骤4.2:计算出每个块的直方图;
步骤4.3:对每个子区域进行灰度直方图“剪切”;
步骤4.4:对每个块在进行直方图均衡化;
步骤4.5:把各个块用线性插值法连接成一幅图像,即为变换后的图像。
进一步的,所述步骤5具体为:
步骤5.1:将血管每15度旋转一次,得到从0度到180度的12个匹配滤波器;
所述的二维匹配滤波器的高斯核是一个二维矩阵,假设q=[x,y]是这个高斯核矩阵中的离散点,φi表示n个匹配滤波器中的第i个高斯核的旋转角度。在计算高斯核旋转矩阵中的系数时,假设旋转中心在(0,0)点处,旋转矩阵为:
旋转后的坐标系中的点为:
步骤5.2:计算匹配滤波器的高斯卷积核;
由于高斯曲线两端沿x轴正负方向无限扩展,需要对其进行截断:u=±3σ,|v|≤L/2。因此,第i个高斯核矩阵中的系数为:
其中,N={[u,v]||u|≤3σ,|v|≤L/2}。
式中,K(x,y)为核函数,u、v分别为在x轴和y轴上的值,σ是高斯函数距离x轴坐标中心的偏移量。σ越大,函数值越向x轴延伸,反之,σ越小,函数值就越向x轴收缩。在y轴方向上被截断的血管长度由L表示。
匹配滤波时,为了使图像原来的背景灰度特性不发生改变,要求卷积核系数平均值为零,因此,匹配滤波器的高斯核的卷积核为:
Ki'(x,y)=Ki(x,y)-mi
其中:
为高斯核的均值,A为N中的点数。
步骤5.3:匹配滤波器设计好之后,利用该滤波器对图像进行卷积运算,然后选择其中最大的卷积值作为增强图像的像素值。假设f(x,y)为原始图像,f'(x,y)为进行匹配滤波增强之后的图像,用n个高斯匹配滤波器对图像中的血管增强的计算公式为:
其中u、v分别为x轴和y轴上的值。
进一步的步骤6所示:
步骤6.1所示:将血管每5度旋转一次,得到从-90度到90度的36个匹配滤波器
所述的二维gabor滤波器的高斯核也是一个二维矩阵,在计算高斯核旋转矩阵中的系数时,假设旋转中心在(0,0)点处,旋转矩阵为:
旋转后的坐标系中的点为:
步骤6.2所示:计算匹配滤波器的高斯卷积核;
步骤5.3匹配滤波器设计好之后,利用该滤波器对图像进行卷积运算,选其中最大的卷积核值作为增强图像的像素值。假设f(x,y)为原始图像,f'(x,y)为进行匹配滤波增强之后的图像,用n个高斯匹配滤波器对图像中的血管增强的计算公式为:
仅使用高斯滤波在微弱血管区域响应上有一定的优势,但是响应产生的背景噪声也会较多,仅使用Gabor滤波器组,微弱血管区域响应不足,且容易产生虚假响应,因此将两者融合,各自发挥自身优势。如图14所示。
进一步的步骤7所示:
将上述分别用高斯滤波之后的结果与经过Gabor滤波之后的结果融合即图像相加,如图15所示,通过对比选用高斯滤波之后的结果与Gabor滤波之后的结果6:4进行融合。
进一步的,所述步骤8具体为:
步骤8.1:将经过三个通道后的图像中的每个像素点减去经过CLAHE与匹配滤波融合处理后的图像对应位置的像素点;
步骤8.2:进行取反运算,即1-F(x,y);
进一步的,所述步骤9具体为:
步骤9.1:使用3×3大小的掩膜进行开运算。
步骤9.2:使用3×3大小的掩膜进行闭运算。
步骤9.3:使用3×3大小的掩膜腐蚀操作。
步骤9.4:将处理后的掩膜乘以处理后的视网膜血管图像得到最终预处理图像11
步骤10:将带有快速连接机制的PCNN与区域生长思想相结合
进一步的,步骤10中所述的PCNN简化模型是20世纪90年代形成和发展的与传统人工神经网络有着根本不同的新型神经网络,它是根据猫的大脑皮层中神经元同步发放脉冲这一现象提出的,在图像分割、边缘检测、降噪、特征提取等领域中有着广泛的应用。PCNN应用于眼底图像处理时,其神经元数目与图像像素数目一致,各神经元与各像素一一对应。PCNN模型可以分为三个部分:输入域、调制域和脉冲产生域。输入域又可以分为:反馈输入域和连接输入域。如图12所示。
进一步的,PCNN用迭代公式表示如下:
Fij[n]=Sij (1)
Uij[n]=Fij[n]{1+βLij[n]} (3)
式中,Sij为外部刺激,即点(i,j)对应像素的灰度值,Fij为神经元的反馈输入项,Lij为神经元的连接输入项,表示一个神经元的周围神经元的输出进行加权求和,Uij为内部活动项,Yij为神经元的脉冲输出,θij为动态阈值。VL为连接输入域的放大系数,β为神经元之间的连接强度系数,Vθ和αθ分别为动态阈值的放大系数和衰减常数因子,Wijkl为连接加权系数矩阵。
通过将内部活动项与动态阈值进行比较,来决定神经元是否输出脉冲,当Uij>θij时,神经元被激发输出脉冲,同时会使得阈值增大,以阻止该神经元再次被激发,不断地迭代,直到某一时刻Uij再次大于θij时,神经元被再次激发,输出脉冲。
基于区域增长的PCNN算法(RG-PCNN)对简化PCNN模型进行了一些改进,在耦合部分如公式(9)连接系数β被集合βn取代:
Uij[n]=Fij[n](1+βn(t)Lij[n]) (6)
在每次迭代中,β值是变化的,先给β设定一个初值,使得第一次发出脉冲的神经元能够就近捕获下一个神经元,然后反复的用设定的较小的δβ来增加它的值。将β变化前后激发的神经元对应的像素区域的边缘比例以及面积比例分别进行比较,如果小于某个阈值,就继续增加β,否则结束迭代过程。
本发明使用单程条件来解决算法的终止,即每个神经元允许发放脉冲,状态从0变1只有一次,算法过程中,将一个区域对应的神经元同时激发。发放脉冲的迭代次数作为区域标记,迭代次数存储在一个矩阵P中,所有元素的初始值设置为0。
为了实现时间独立性,本发明采用发出脉冲的神经元确定的阈值Wn,n表示迭代次数。其他神经元的阈值设置一个很大的阈值Ω,以防止同一位置重复激发。
将Wn设置为没有激发状态神经元对应的像素的最大值,即将一个强度最大的没有激发过的神经元的初始状态设置为激发状态。如果没激发的神经元具有相同的最大强度值,则只选择一个激发。
由于在PCNN模型中连接域的脉冲传输具有延时性,这样会导致图像中出现断点或者产生不连续的区域。为了减少这些影响,采用快速连接模型,它可以使得已经激发的神经元所引起的自动波在一次迭代过程中充分传播,进而使所有被捕获的神经元可以在同一次迭代过程中被激发从而产生脉冲。
为了得到更好的血管分割结果,本文选取以下几种情况作为局部迭代及全局迭代终止条件:
(1)所有的神经元都被激发;
(2)β大于βmax;
(3)生成血管边缘像素点数量占血管像素点总数的比例e小于设定的阈值;
(4)血管占整幅图的面积比m大于设定阈值;
(5)快速连接终止。
其中(1)为全局迭代条件,(2)(3)(4)(5)为局部迭代条件。
条件(1)即整幅图片中每个像素都被划分到不同区域,不存在没有划分的像素,该条件可以通过判断P中有没有为0的元素来确定。因为P中的初始值都为0,每经过一次迭代,迭代次数赋给激发神经元对应的像素;条件(2)即每次迭代的连接系数强度β大于给定值βmax则终止。相比于原始RG-PCNN,本文根据实际应用场景增加了(3)(4)两个终止条件。其中(3)是为了判断当前种子生长合适的生长区域,如果种子生长的血管边沿像素点总个数与血管总像素点个数的比值小于设定的阈值e则删除掉这部分伪血管,终止本次迭代,转而从剩下未激发的像素点中选择最大的种子点继续上述步骤。(4)是限制血管生长总面积占整幅图的比例小于设定的阈值m,防止血管过生长。条件(5)即没有神经元再被捕获。比较迭代前后激发的神经元个数,如果相等,则停止,如果迭代后激发的神经元数量仍在增加,则迭代继续。
进一步的,计算等式(1)-(5)式进行第一次迭代,重复计算:
Uij[n]=Fij[n]{1+β(n)Lij[n]} (10)
直到Y不再发生变化为止,然后按照公式(12)更新阈值,开始新一轮的迭代:
所述步骤10具体为:
步骤10.1:设置连接强度系数初始值β、连接系数的增加项δβ、连接强度系数最大值的初始值βmax,生成血管边缘像素点数量占血管像素点总数的比例e,血管占整幅图的面积比m;停止阈值T=255;
步骤10.2:将预处理结果进行归一化Sy,作为PCNN的外部输入Fij[n]=Sy;
步骤10.3:选取的种子点,将Sij中最大的像素点作为种子点将它的值赋为1,Y(ijseed)=1;
步骤10.4:设定PCNN阈值Wn取没有激发状态神经元对应像素的最大值,即Tij[n]=Wn,其他神经元阈值被设置为一个很大的常数Ω(默认为500);
步骤10.5:引入快速连接机制进行迭代分割,即定义一个与PCNN_Y同阶的矩阵Y0,使得Y0=1,当β≤βmax时,如果Y0与PCNN_Y中存在不相等的值,则Y0=PCNN_Y,然后引入快速连接机制进行迭代,重新计算PCNN_Y;
步骤10.6:反复循环,不断地比较迭代后的PCNN_Y和迭代前的值Y0,直到两者相等为止;
步骤10.7:用较小的常数δβ来增加它的值,并将Y0的值再赋为1,开始新的循环;
步骤10.8:重复上述过程,直到β>βmax为止,血管完成第一轮生长;
步骤10.9:最后通过计算生成血管边缘像素点数量占血管像素点总数的比例e小于等于0.2时并且血管占整幅图的面积比m小于等于0.11时,生长无效,删除掉一部分伪血管,开始新一轮地生长,并设置了T=255作为停止条件。
整个算法描述如下:
设置β0、δβ、βmax、T、e、m;
将预处理结果归一化为Sy,作为PCNN的外部输入:Fij[n]=Sy;
选取的种子点,将Sij像素值中最大的点对应像素值赋为1,即Y(ijseed)=1;
找到没有激发神经元对应像素值中最大值置为阈值Wn,即Tij=Wn,其他神经元的阈值被设置为一个很大的常数Ω(默认为500)
whileβ≤βmax
while PCNN_Y!=Y0
Y0=PCNN_Y
Uij[n]=Fij[n]{1+β(n)Lij[n]}
end。
β=β+δβ
Y0=1
end
统计生成血管边缘像素点比例e,血管占整幅图的面积比m:
if sum(e>0.2&&m>0.11)
n=n+1
统计P==0索引值
end
本发明的有益效果是:
(1):从DRIVE数据库中随机抽出10幅图,通过准确度(Acc),灵敏度(Sen),特异性(Spe)的计算,反映出预处理中使用高斯滤波与Gabor滤波相结合的办法,比仅使用Gabor滤波或高斯滤波处理结果好。
表1 Gabor滤波的结果
表2高斯滤波的结果
表3 Gabor滤波与高斯滤波融合的结果
(2):提供了一种基于区域增长PCNN的视网膜血管分割方法。本发明在对视网膜血管预处理的基础上首先选择出一定的血管区域作为初始种子区域;接着,将带有快速连接机制的PCNN与种子区域增长思想相结合,通过设置PCNN中的连接强度系数最大值和停止条件,如果生长无效则会删除掉那部分伪血管,解决了以前方法中断点不能继续生长的问题,实现眼底图像中血管区域的自动生长,从而有效地提取出视网膜图像中的血管。