一种改进NAS-RIF的OCT图像盲复原方法与流程

本发明涉及一种光学相干层析(opticalcoherencetomography,oct)图像盲复原方法，属于视觉检测领域，尤其涉及一种改进非负性和递归逆滤波(nonnegativityandsupportrecursiveinversefiltering,nas-rif)的oct图像盲复原方法。

背景技术

oct技术利用光源低相干性获取待测样品内部结构信息，已应用于生物医学、临床医学和物体无损检测等领域，然而oct系统在成像过程中易受到噪声的影响，造成图像退化，掩盖和模糊图像的细节信息，直接影响着oct图像分析精度。因此，实现oct图像的复原，获得结构清晰、边缘明显的高信噪比oct图像，能够提高oct图像分析与检测的准确性。

现有的oct图像复原方法主要有非盲复原法和盲复原法。非盲复原法根据实际oct图像的点扩散函数(pointspreadfunction,psf)实现oct图像复原，然而实际应用中oct图像的psf较难准确获得，导致采用非盲复原法难以获得高质量的oct复原图像；而盲复原法中的nas-rif方法能在psf未知的情况下，利用oct图像的目标支持域和非负性，施加限制条件约束复原过程，实现可靠的图像复原，然而该方法在oct图像复原的过程中易受观测图像噪声的影响，使得估计图像产生偏差，噪声鲁棒性低，制约了oct图像复原质量的提高。因此，通过改进nas-rif方法，提高oct图像盲复原过程中的噪声鲁棒性，实现高质量的oct图像复原。

技术实现要素：

本发明的目的是提高nas-rif方法在oct图像盲复原过程中的噪声鲁棒性，获得高质量的oct复原图像。首先，应用几何非线性扩散(geometricnonlineardiffusionfilter,gnldf)算法对oct图像预处理，并采用基于泊松分布的最小误差阈值估计图像的目标支持域；然后，通过引入oct估计图像的梯度信息，构建基于各向异性总变分(anisotropictotalvariation,atv)正则化的改进代价函数，并采用分裂bregman迭代算法优化改进代价函数，进而调整fir滤波器系数，使投影图像不断接近真实图像，实现oct图像的盲复原。

本发明采用的技术方案为一种改进nas-rif的oct图像盲复原方法，该方法具体包括以下步骤：

步骤一：对oct原图像进行gnldf滤波，得到gnldf滤波后的图像；

步骤二：将gnldf滤波后的图像经过fir滤波器，得到估计图像然后采用基于泊松分布的最小误差阈值估计图像的支持域，再经过nl滤波器，得到投影图像

步骤三：构建基于atv正则化的改进代价函数j，并采用分裂bregman迭代方法对改进代价函数进行优化；

步骤四：重复步骤一到步骤三，直至满足终止条件时停止迭代，此时的投影图像即为盲复原后的oct图像。

所述步骤一，具体包括：

oct图像为i，利用式(1)估计出参数δ。

式中，ii,j代表坐标为(i,j)处的像素值；mean(·)为图像像素的均值；median(·)为中位数。

利用估计出的参数δ和待处理像素点is，以及is右边和左边相邻的像素点的强度ie和iw，由式(2)得到参数dx，由式(3)得到x方向像素平均值ax，在此基础上，由式(4)得到参数i's,x，从而得到参数px。

式中，|·|为绝对值。

ax＝(ie+iw)/2(3)

px＝i's,x-ax(5)

利用参数dx和px，由式(6)得到x方向的扩散系数c(dx,px)。

y方向的扩散系数c(dy,py)，进而由式(7)得到离散形式的gnldf。

式中，is^t是在t时刻的扩散图像，即为滤波后的图像g(x,y)；δt为时间步长大小；p＝e,w,n和s；e、w、n和s分别代表了中心像素点的右、左、下和上像素点。

所述步骤二，具体包括：

根据式(8)，将滤波后的图像g(x,y)经过fir滤波器u(x,y)，得到估计图像

式中，*为卷积。

根据式(9)和式(10)分别估计图像背景区域和支持域的先验概率p0(t)和p1(t)，平均值μ0(t)和μ1(t)。

式中，imax为最大的像素值；h(ii)为图像直方图归一化处理结果；ii表示像素值大小。

在此基础上，根据式(11)得到最佳阈值t。

根据获得的最佳阈值t，将估计图像的像素分为背景和目标区域两类，大于t的像素构成目标支持域，小于t的像素构成背景区，组成与退化图像尺寸相同的二维数组模板b(x,y)。

式中，1表示支持域；0表示非支持域。

根据式(12)得到的b(x,y)，将估计图像经过nl滤波器，对图像非负性进行约束，支持域内大于0的像素大小保持不变，小于0的置于0，支持域外的像素大小置为背景平均值lb，得到投影图像

式中，为投影图像在坐标(x,y)处的像素值大小；dsup为支持域内像素点的集合；为支持域外像素点的集合；lb为图像背景灰度均值。

所述步骤三，具体包括：

利用由步骤二获得的估计图像和投影图像根据式(14)和式(15)得到改进代价函数的第一部分和第二部分j1、j2。

j2＝γ(∑u-1)²(15)

式中，μ为平衡参数；常数γ只有在lb为零值时不取零，其余均为正数。

然后根据式(16)和式(17)，计算出估计图像x方向和y方向的梯度和

式中，m×n为的大小；代表坐标(x,y)处的像素强度。

根据估计图像的梯度，进一步由式(18)得到改进代价函数的第三部分jatv。

式中，||·||1代表l1范数。

利用获得的j1、j2和jatv，根据式(19)得到基于atv正则化的改进代价函数j。

j＝jatv+j1+j2(19)

根据式(8)，将改进代价函数改写为

采用分裂bregman迭代方法对改进代价函数(20)进行优化，将式(20)转变成无约束最优化问题和bregman迭代，得到迭代后的u(x,y)，具体迭代过程如下：

引入辅助变量dx和dy，令其分别为式(21)和式(22)

则式(20)转化为式(23)所表示的约束问题

通过借助辅助变量dx、dy和二次惩罚项，将式(23)转化为

式中，λ>0，为惩罚因子；控制惩罚项的权重；||·||2代表l2范数；bx和by为引入的两个变量。

再利用迭代次数为k时的u^k、和通过bregman迭代，由式(25)得到迭代次数为k+1时式(24)中的u^k+1、和

式中，shrink(：，：)为阈值收缩算子，定义为

shrink(ui,δ)＝sign(ui)max(0,|ui|-δ)(26)

分裂bregman迭代通过引入辅助变量、变量分离和bregman迭代步骤，将具有l1项函数的优化问题式(20)转化为一系列无约束最优化问题式(25)，从而实现改进代价函数的优化。

所述步骤四，具体包括：

重复步骤一到步骤三，直至满足式(27)或式(28)时停止迭代，此时的投影图像即为盲复原后的图像。

k≥kmax(27)

式中，k为迭代次数；kmax为设定的最大迭代次数。

j≤ε(28)

式中，ε为设定值，为较小的正数。

本发明的优点为：通过采用gnldf算法对oct图像预处理，有效滤除部分噪声的同时保留了图像的细节信息，减小了噪声对复原过程的影响；通过构建的基于atv正则化的改进代价函数，充分利用了oct估计图像的梯度信息，提高了nas-rif方法在oct图像复原过程的噪声鲁棒性，且能够保留图像的细节信息，获得了结构清晰、边缘明显的高信噪比oct复原图像。

附图说明

图1是本发明所述的一种改进nas-rif的oct图像盲复原方法的构建框图。

图2是具体实施例中的扫频oct系统对橘子组织的成像结果以及处理结果图，其中，(a)为采用扫频oct系统对橘子组织进行扫描成像结果；(b)是采用gnldf方法对(a)去噪后的结果；(c)是采用基于泊松分布的最小误差阈值法支持域提取结果；(d)是本发明方法对(a)复原后的结果。

图3为图2中各图的局部放大图。

具体实施方式

下面结合实例及附图对本发明作进一步的描述，需要说明的是，实施例并不限定本发明要求保护的范围。

实施例

采用中心波长为1310nm，带宽为170nm，快扫频率为50hz，轴向分辨率为10μm，横向分辨率为20μm的扫频oct设备对橘子组织进行扫描成像，获得的图像如图2(a)所示。

将本发明应用到上述橘子组织的oct图像中，具体步骤如下：

步骤一：利用式(1)至式(6)，计算出oct图像x方向的扩散系数c(dx,px)和y方向的扩散系数c(dy,py)，利用式(7)计算出t时刻的扩散图像is^t，其中δt为0.25，t为10，此时的扩散图像is^t即为滤波后的图像g(x,y)。

步骤二：将滤波后的图像g(x,y)经过5×5大小的fir滤波器u(x,y)，根据式(8)得到估计图像利用式(9)和式(10)得到估计图像背景区域和支持域的先验概率p0(t)和p1(t)，平均值μ0(t)和μ1(t)，再利用式(11)得到最佳阈值t；

将估计图像的像素分为背景和目标区域两类，大于t的像素构成目标支持域，小于t的像素构成背景区，根据式(12)组成与退化图像尺寸相同的二维数组模板b(x,y)；

将估计图像经过nl滤波器，根据得到的b(x,y)，利用式(13)对图像非负性进行约束，支持域内大于0的像素大小保持不变，小于0的置于0，支持域外的像素大小置为背景平均值lb，得到投影图像

步骤三：利用式(14)和式(15)，由估计图像和投影图像得到改进代价函数的第一部分和第二部分j1、j2，其中μ为200，γ为10，利用式(16)和式(17)计算出估计图像x方向和y方向的梯度和从而根据式(18)得到改进代价函数的第三部分jatv，根据得到的j1、j2和jatv，利用式(19)构建基于atv正则化的改进代价函数j；

利用式(20)到式(24)引入辅助变量dx和dy，并进行变量分离，根据式(25)对改进代价函数j进行优化，更新fir滤波器系数，其中λ为10。

步骤四：重复步骤一到步骤三，直到迭代次数大于10次，停止迭代，此时的投影图像即为盲复原后的oct图像。

图像质量评价指标为信噪比snr，噪声对比度cnr和等效视数enl，cnr可以衡量目标区域和背景噪声区域的对比度，enl衡量受噪声污染区域的图像平滑度，snr和cnr越大表示图像的去噪效果越好，enl越大表示处理后的图像边缘越平滑，评价指标对比结果如表1所示。

在图中选择r个目标区域以及一个背景区域进行计算，其中目标区域用白色点划线框标出，背景区域用白色实线框标出，如图2(a)所示，本发明中选取3个目标区域。式中i代表图像的像素值，σ为背景区域的标准差，μr和σr分别是第r个目标区域的均值和标准差，μb和σb分别是背景区域的均值和标准差。

表1oct图像复原结果评价指标对比

图2(a)和图3(a)为oct原图，从图中可以看出噪声严重且图像边缘不清晰；图2(b)和图3(b)为经过gnldf方法预处理后的图像，从图中可以看出部分噪声得到了有效地滤除，且图像边缘保持效果较好；图2(c)和图3(c)为采用基于泊松分布的最小误差阈值法提取的支持域，从图中可以看出没有损失图像的细节信息，目标区域提取结果准确；图2(d)和图3(d)为采用本发明方法处理后的图像，从图中可以明显看出，噪声得到了更加有效的抑制，图像边缘更加清晰，由于在代价函数中引入了抑制噪声的正则化项，因此复原的图像边缘相对平滑，同时从表1中可看出，snr、cnr和enl相较于原图均有了明显提高，进一步表明了本发明方法有效地提高了nas-rif方法在oct图像复原过程的噪声鲁棒性，且能够保留图像的细节信息，获得了结构清晰、边缘明显的高信噪比oct复原图像。