一种自适应尺度分解的眼底图像血管随机共振检测方法与流程

本发明属于数字图像处理领域，主要涉及一种自适应尺度分解的眼底图像血管随机共振检测方法。

背景技术：

眼底图像血管的自动检测，将有助于改善眼部疾病临床辅助诊断的效率和准确率。传统分割方法将噪声视为干扰信号，试图通过抑制噪声的手段来实现图像的分割，但在噪声处理的过程中，不可避免的会影响到低对比度弱信号的检测，因此它们不适合类似眼底图像血管自动检测的应用场景。目前随机共振方法成为了弱信号检测的研究热点，它认为噪声的存在对弱信号的检测是有利的，因为随机共振方法可以实现噪声能量向有用弱信号能量的转移，从而能够改善眼底视网膜微弱边缘血管的检测。但目前的随机共振方法，大多是对整幅图像进行一种空间上的操作，而由于眼底图像中的血管呈现出多对比度分布特性，因此所获得的结果只能是一种全局意义上的最优随机共振响应。

技术实现要素：

本发明利用非下采样轮廓波变换(nsct)对源图像进行尺度分解，分解为高频子带图像和低通子带图像，其中低通子带图像主要反映了源图像的整体轮廓，高频子带图像则携带了源图像的细节信息。具体对于眼底视网膜血管图像，低通子带图像主要体现高对比度血管信号的特性，而高频子带图像中主要包含低对比度微细血管信号和噪声信号。为了实现低对比度微细血管的有效检测，本发明对分解得到的高频子带图像进行最优随机共振调制，并通过非下采样轮廓波变换分解尺度的自适应控制，从而改善对微细血管的检测性能。最后对不同尺度下获得的低通子带图像以及最优随机共振响应进行融合，并通过非下采样轮廓波逆变换重构得到眼底图像血管的检测结果。

本发明包括以下步骤：

步骤1：对眼底图像进行第1级nsct变换，得到1个低通子带图像和8个方向的高频分量子带图像。

步骤2：基于能量最大规则，对步骤1中得到的8个方向的高频分量子带图像进行融合，将融合后的子带合成图作为duffing随机共振模型的输入信号，并且以峰值信噪比为评价指标调节duffing随机共振模型中的附加内噪声强度，得到峰值信噪比下的最优随机共振响应。

步骤3：根据步骤2中最优随机共振响应时所对应的峰值信噪比指标，若其高于设定的最低标准值t，则对步骤1得到的低通子带图像继续进行下一级的nsct变换，此时将得到新尺度下的1个低通子带图像和8个方向的高频分量子带图像。

步骤4：重复执行步骤2和3，直到峰值信噪比指标低于设定的最低标准值t，假设此时进行了第n级nsct变换，则结果包括不同尺度下的n个低通子带图像以及n个最优随机共振响应。

步骤5：基于区域能量加权平均规则，对步骤4得到的n个低通子带图像进行融合；同时基于峰值信噪比数值的加权方法，对步骤4得到的n个最优随机共振响应进行加权融合。

步骤6：对步骤5得到的融合结果进行nsct逆变换重构，并通过二值化阈值处理后得到最终的眼底图像血管检测结果。

本发明具有的特点为：本发明对眼底图像进行非下采样轮廓波变换，由于分解得到的高频子带图像主要反映了眼底图像低对比度血管分布的边缘信息，因此对高频子带图像进行针对性的最优随机共振调制，改变了传统随机共振方法对整幅空域图像进行全局优化的伪最优思路，从而能更好地检测出眼底图像中的低对比度微细血管信号。本发明根据高频子带图像最优随机共振响应所对应的峰值信噪比指标，来自适应控制非下采样轮廓波变换的分解尺度；分别基于区域能量加权平均规则和基于峰值信噪比数值的加权方法，对不同尺度下获得的低通子带图像以及最优随机共振响应进行融合；最后对融合结果进行非下采样轮廓波逆变换重构，并利用二值化阈值处理得到血管检测结果图。本发明所提方法适用于眼底图像多对比度血管的检测场景。

附图说明

图1为三级分解滤波器组示意图；

图2为nsdfb八方向划分示意图；

图3为nsct结构分解示意图；

图4为本发明专利的方法流程示意图。

具体实施方案

本发明一种自适应尺度分解的眼底图像血管随机共振检测方法，该方法具体包括以下步骤：

步骤1：将原始眼底图像记为源图像p，m和n表示源图像p的行列尺寸。对p进行第1级nsct变换，得到1个低通子带图像和8个方向的高频分量子带图像。nsct具有平移不变特性使得这些子带图像与源图像具有相同的尺寸大小m×n。

本发明中的nsct变换首先进行非下采样金字塔分解(nsp)，通过双正交小波‘9-7’构建滤波器组得到一个低通滤波器l和高通滤波器h，设则滤波器l和h的矩阵构造分别由滤波器系数li,j和hi,j决定，取c＝3，滤波器l和h的矩阵形式分别如式(1)和(2)所示。

其中l0,0＝0.3245066567828517，l0,1＝l1,0＝0.1643662782633664，l0,2＝

l2,0＝-0.18830187186854155，l0,3＝l3,0＝-0.013629023704276572，

l1,1＝0.0740121520471188，l1,2＝l2,1＝-0.02131840020457929，

l2,2＝-0.01306636838062783，l1,3＝l3,1＝-0.007296122436400075，l2,3＝

l3,2＝0.001762583572331734，l3,3＝0.003080092825681339，

h0,0＝1.669154993235563，h0,1＝h1,0＝-0.24975066285895303，h0,2＝

h2,0＝-0.055805733503531366，h0,3＝h3,0＝0.028522425973203154，

h1,1＝-0.1692097602584649，

h1,2＝h2,1＝-0.022942609349085995，h2,2＝0.01351698519464063，h1,3＝

h3,1＝0.021391819479902367，h2,3＝h3,2＝0.008556727791960947，

h3,3＝0.0014261212986601578。

源图像p分别与低通滤波器l和高通滤波器h进行卷积运算，如式(3)和(4)所示，得到一个低通子带图像e1和一个高频分量图像f1。其中*为卷积符号，下同。

e1＝p*l(3)

f1＝p*h(4)

获得的高频分量图像f1可通过非下采样方向滤波器组(nsdfb)进行多方向分解，如式(5)所示。其中三级分解滤波器组进行方向分解的具体流程示意图如附图1所示，8个方向的选择示意图如附图2所示。采用梅花采样方向滤波器组，加入扇形垂直滤波器a和扇形水平滤波器b，对高频分量图像进行三级分解，可以得到8个方向的高频分量子带图像。其中·为矩阵相乘符号。

2级分解梅花采样矩阵定义为：

3级分解梅花采样矩阵定义为：

扇形滤波器矩阵定义为：

即对源图像p进行一次nsct变换后得到1个低通子带图像e1和8个方向的高频分量子带图像f1,1～f1,8。

步骤2：基于能量最大规则，按式(9)对前述步骤得到的8个方向的高频分量子带图像进行融合，得到融合后的子带合成图f1max。

f1max(x,y)＝max(f1,j(x,y),j＝1,2,...8)(9)

其中f1,j(x，y)表示高频分量子带图像f1,j中(x，y)位置所对应的像素值，其中0<x≤m，0<y≤n，max为选取最大值的函数。

利用逐行扫描将f1max降维成一维信号，然后将其作为duffing随机共振模型的输入信号，duffing随机共振模型可以用式(10)所示的方程描述：

其中，k为阻尼比惯性系数；s(t)表示随机共振模型的输入，即为f1max降维所对应的一维信号；ξ(t)为附加内噪声，且满足：<ξ(t)>＝0，<ξ(t)ξ(t')>＝2dδ(t-t')，d为噪声强度；v表示随机共振模型的输出信号。u(v)为如式(11)所示的势函数，其中参数a和b决定了系统的势阱位置和势垒高度。参数设定如下：k＝1，a＝1，b＝1。

以峰值信噪比为评价指标调节duffing随机共振模型中的附加内噪声强度，将输出信号v逐行重构回二维信号得到随机共振响应图像y1，当y1的峰值信噪比达到最大值时，表明此时的附加内噪声强度对应于最优随机共振响应。其中峰值信噪比psnr的计算方式如式(12)和(13)所示：

步骤3：根据步骤2中最优随机共振响应时所对应的峰值信噪比指标，若指标高于最低标准值t，则将前述步骤中得到的低通子带图e1继续进行下一级的nsct变换，此时将得到新尺度下的1个低通子带图像e2和8个方向的高频分量子带图像f2,1～f2,8。其中最低标准值t的设定，可根据源图像p的局部区域信噪比结果进行自适应调整，一般可随机选取源图像p的若干个局部区域并计算它们所对应的信噪比，比较获得最小的信噪比数值，最后将t设定为此数值的50％。

步骤4：重复执行步骤2和步骤3直到峰值信噪比指标低于设定的最低标准值t，则计算自动结束。假设此时进行了第n级nsct变换，则结果包括不同尺度下的n个低通子带图像e1,e2,...en以及n个最优随机共振响应结果y1,y2,...yn。

步骤5：对步骤4中得到的n个不同尺度的低通子带图像e1,e2,...en通过区域能量加权平均规则进行融合。图像en的区域能量enn通过式(14)和(15)定义。

其中en(i,j)为en在(i，j)点的像素值，enn(i,j)为en在(i，j)点处的能量，wdα,β为对应的窗口函数。考虑到低通子带图像的方向性和中心像素特点，选定α＝{-1,0,1}，β＝{-1,0,1}，设定窗口函数wdα,β如式(16)所示。

计算各个低通子带图像e1,e2,...en的区域能量en1,en2,...enn，则低通子带图像ei的权值为：

n个低通子带图像e1,e2,...en经过式(18)进行融合，得到低通子带融合图像e^*。

基于峰值信噪比数值的加权方法对步骤4中得到的n个最优随机共振响应y1,y2,...yn进行加权融合，则yi的权值为：

其中psnri，i＝1,…n，表示第i次最优随机共振响应的峰值信噪比数值。

n个最优随机共振响应y1,y2,...yn经过式(20)进行融合，得到融合结果y^*。

步骤6：对步骤5中得到的融合结果e^*和y^*进行nsct逆变换重构，得到重构后的图像，图像经过二值化阈值处理后即可得到最终的眼底图像血管检测结果。上述步骤设计中所用到的nsct结构分解示意图如附图3所示。本发明专利的流程示意图如附图4所示。