技术领域
[0001] 本
发明涉及一种图像去噪方法，特别是涉及一种阈值寻优的高保真各向异性滤波方法，属于
图像处理技术领域。
背景技术
[0002] 图像去噪是图像处理与计算机领域中最基本的问题之一。图像去噪是图像的
边缘检测、增强和分割等图像处理中的一个非常重要预处理。现阶段，图像去噪方法种类繁多，有基于偏微分方程、概率论、小波等图像去噪方法。
[0003] 近20年来，基于偏微分方程的方法在图像去噪中得到了广泛的应用。Perona和Malik首次提出了经典的
各向异性扩散模型，即PM模型，该模型由一个关于图像梯度模值的扩散函数控制扩散强度；自PM模型以后，基于PDE的去噪方法继续发展，其有效地提高了边缘检测能
力。
[0004] 目前，Zhang K B，Gao X B等人提出了一种NLM
滤波器，该方法较好地保持了图像纹理和细节特征，但是计算复杂高，处理速度慢。BM3D去噪方法主要通过
块匹配以及三维变换域滤波等技术进行两次去噪，第一次
基础去噪为第二次最终去噪提供权值，使得最终去噪效果更加好，该方法去噪后的图像不仅有较高的
信噪比，视觉效果也很好。Deledalle等提出了PPB滤波器，该滤波器通过重新定义
像素之间的相似度，并用加权最大似然估计完成图像的重建，在
合成孔径雷达(SAR)图像的处理上达到较好的效果，但实验结果表明该滤波器运行时间较长，而非
迭代滤波方法不能保留图像纹理等细节特征，SAR图像噪声的先验信息未充分利用。
[0005] 小波具有良好的时频变换特性，使得它在图像去噪领域发挥着重要作用，在基于
小波变换的去噪方法中，小波硬阈值去噪
算法和小波软阈值去噪算法是较为有效地降噪手段。目前，基于小波的图像阈值去噪取得了很多研究成果，在扩散滤波以及它的相应变化方法方面，人们逐渐对
离散小波变换和偏微分方程之间的联系产生了兴趣。
[0006] 在各向异性扩散滤波方法中，边缘检测的准确性对滤波结果有重要的影响。然而，目前基于偏微分方程去噪方法的扩散强度大都使用梯度信息检测边缘，在图像的尖峰、窄边缘处一阶偏导数为0，对图像边缘纹理的尖峰进行了平滑处理，磨平了图像的尖峰、
角点等细节信息，从而使图像失真。
[0007] 总之，上述传统方法时效性低，复杂度较高，且在处理过程中图像降噪不稳定，有明显的“阶梯”效应，并且保边缘性效果不是很好。
发明内容
[0008] 本发明的主要目的在于，克服
现有技术中的不足，提供一种阈值寻优的高保真各向异性滤波方法，可有效保护图像的尖峰、角点等细节信息，提高边缘检测的准确性，增强图像滤波效果，去噪性能优异；而且运行时间短，有利于实际应用。
[0009] 为了达到上述目的，本发明所采用的技术方案是：
[0010] 一种阈值寻优的高保真各向异性滤波方法，包括以下步骤：
[0011] 1)输入图像，采用各向异性扩散模型、即PM模型对图像进行预处理后输出含噪声的噪声图像；
[0012] 2)判断是否考虑局部特征处理，若是、则进入步骤3)，若否、则进入步骤6)；
[0013] 3)用小波变换提取噪声图像的高频部分，在高频部分用二阶微分量
曲率模值m来反映局部信息，并建立高保真各向异性滤波模型；
[0014] 根据式(1)计算曲率模值m，
[0015]
[0016] 式(1)中，Ixx、Ixy分别为Ix在x,y方向的方向导数，Iyy为Iy在y方向的方向导数，其中，Ix为图像I在x方向分量，Iy为图像I在y方向分量；
[0017] 高保真各向异性滤波模型为式(2)，
[0018]
[0019] 式(2)中，x,y分别为图像的行数与列数，t为时间，div为散度算子，W为小波变换，×为用小波对图像做小波分解以提取图像高频部分，▽为梯度算子，m为曲率模值，g(·)为依赖于图像的单调递减的扩散函数，I(x,y,t)为t时刻第x行第y列的像素值、构成平滑图像， 为像素值I(x,y,t)对时间t的导数，I(x,y,0)为原始图像，I0(x,y)为0时刻原始图像的像素值；
[0020] 其中，扩散函数表达式采用 x为检测算子，k为阈值；
[0021] 4)用最小均方算法进行阈值寻优，进一步控制扩散函数的扩散强度；
[0022] 根据最小均方算法中的
梯度下降法来求取最优阈值k，表达式为式(3)，[0023] k(n+1)＝k(n)-μ▽k(n) (3)
[0024] 式(3)中，k(n+1)为下一时刻的阈值，k(n)为当前时刻的阈值，μ为步长，▽k(n)为均方误差的梯度；
[0025] 5)用建立的高保真各向异性滤波模型对提取的高频部分进行处理，对处理后的高频系数和原来的低频系数进行小波重构，得到去噪后图像并输出；
[0026] 将平行于y方向的边缘，忽略y方向的影响和小波变换且Ix＞0，则反映局部信息的检测算子简化为 其中，u＝|W×(▽I+m2)|为检测算子；
[0027] 此时高保真各向异性滤波模型的简化模型为式(16)，
[0028]
[0029] 式(16)中，
[0030]
[0031] 6)输出图像。
[0032] 本发明进一步设置为：所述步骤4)根据最小均方算法中的梯度下降法来求取最优阈值k，具体为，
[0033] 4-1)对噪声图像进行J层离散小波变换，得到原始
信号的小波系数，记为向量Y，且Y＝[AJ,HJ,VJ,DJ,HJ-1,VJ-1,DJ-1,…,H1,V1,D1]，其中，J为层数，向量A为低频系数，向量H为
水平高频系数，向量D为对角高频系数；
[0034] 设g(Y)是Y的函数，表达式为式(4)，
[0035]
[0036] 式(4)中， 为基于向量Y对信号w的估计值，故g(Y)为RN到RN的映射，RN为N维向量空间；
[0037] 4-2)基于Sure无偏估计进行阈值的选择，通过
风险函数R(t)式(5)来定义，[0038]
[0039] 式(5)中，f、分别为原始信号和原始信号的估计，N为维数；
[0040] 由于小波变换具有
正交性，故风险函数可同样在小波域中表达成式(6)，[0041]
[0042] 式(6)中，X、 分别为原始信号的小波系数和原始信号小波系数的估计；
[0043] 4-3)现记Y的风险函数为式(7)，
[0044]
[0045] 则获得其数学期望为式(8)，
[0046]
[0047] 式(8)中，σn为标准差，V为向量Y与原始信号的小波系数X的差值；
[0048] 4-4)可知，V＝Y-X，当V服从高斯分布时，如下式(9)成立，
[0049]
[0050] 式(9)中，Yi为各尺度下的小波系数向量，Vi为各尺度下的Y向量与X向量的差值，i＝1,2,…,N；
[0051] 则得到式(10)，
[0052]
[0053] 式(9)中，ti为t在各尺度下的值；
[0054] 4-5)将式(10)代入式(8)，得到风险函数的数学期望ER(t)，为式(11)，[0055]
[0056] 根据Sure无偏估计式(11)，可得均方误差 为式(12)，
[0057]
[0058] 式(12)中，为阈值函数作用于Y后得到的信号估计值，gi为各尺度下 与Y的差值；
[0059] 则均方误差 的梯度表达式为式(13)，
[0060]
[0061] 因为 所以得gi为式(14)，
[0062] gi＝X(kn,Yi)-Yi (14)
[0063] 4-6)将式(14)代入式(13)，得式(15)
[0064]
[0065] 根据式(15)计算出▽k(n)，再根据式(3)计算出k(n+1)，从而获得最优阈值。
[0066] 本发明进一步设置为：所述步骤5)中的式(17)，当在噪声图像的拐点时，由Marr边缘
位置知，噪声图像的二阶微分量Ixx＝0，此时有
[0067]
[0068] 又因为 求导得 则得到
[0069]
[0070] 将式(19)代入式(18)中，则得到式(20)，
[0071]
[0072] 式(19)、式(20)中，k为步骤4)寻优到的最优阈值。
[0073] 本发明进一步设置为：所述步骤5)中的式(17)，当处于噪声图像的尖峰、角点时，Ix＝0，此时，反映局部信息的检测算子 则式(17)简化为式(21)，
[0074]
[0075] 式(21)中，k为步骤4)寻优到的最优阈值。
[0076] 本发明进一步设置为：所述步骤6)输出图像之前进行所有像素都直通的判断，若是、则直接输出图像，若否、则返回步骤1)。
[0077] 与现有技术相比，本发明具有的有益效果是：
[0078] 1、在去噪性能方面，利用二阶微分量曲率模值反映图像局部信息，可避免将图像的尖峰、角点等误认为是噪声，有效地保护了图像的尖峰、角点等细节信息，用最小均方算法设计阈值，提高了边缘检测的准确性，增强了图像滤波的效果，去噪性能更具优越性。
[0079] 2、在运行时间方面，通过用二阶微分量曲率模值进行检测算子，并用最小均方算法进行阈值寻优，能快速地进行各向异性滤波，从而减少运行时间，而且又采用AOS算法进行数值分解，并用Thomas方法进行求解，可进一步减少运行时间。
[0080] 上述内容仅是本发明技术方案的概述，为了更清楚的了解本发明的技术手段，下面结合
附图对本发明作进一步的描述。
附图说明
[0081] 图1为本发明阈值寻优的高保真各向异性滤波方法的
流程图；
[0082] 图2为本发明仿真实验所用的原始图像，依次为Nuist图像和Lena图像；
[0083] 图3为原始图像分别采用不同模型的滤波结果(方差σ2＝20)；
[0084] 图4为原始图像分别采用不同模型滤波的边缘提取(方差σ2＝20)；
[0085] 图5为原始图像分别采用不同模型的滤波结果(方差σ2＝50)；
[0086] 图6为原始图像分别采用不同模型滤波的边缘提取(方差σ2＝50)；
[0087] 图7为原始图像分别采用不同模型滤波的运行时间；
[0088] 图8为原始图像分别采用不同模型在不同方差噪声下
峰值信噪比。
具体实施方式
[0089] 下面结合
说明书附图，对本发明作进一步的说明。
[0090] 如图1所示，本发明提供一种阈值寻优的高保真各向异性滤波方法，包括以下步骤：
[0091] 1)输入图像，采用各向异性扩散模型、即PM模型对图像进行预处理后输出含噪声的噪声图像。
[0092] 2)判断是否考虑局部特征处理，若是、则进入步骤3)，若否、则进入步骤6)。
[0093] 3)用小波变换提取噪声图像的高频部分，在高频部分用二阶微分量曲率模值m来反映局部信息，并建立高保真各向异性滤波模型；
[0094] 根据式(1)计算曲率模值m，
[0095]
[0096] 式(1)中，Ixx、Ixy分别为Ix在x,y方向的方向导数，Iyy为Iy在y方向的方向导数，其中，Ix为图像I在x方向分量，Iy为图像I在y方向分量；
[0097] 高保真各向异性滤波模型为式(2)，
[0098]
[0099] 式(2)中，x,y分别为图像的行数与列数，t为时间，div为散度算子，W为小波变换，×为用小波对图像做小波分解以提取图像高频部分，▽为梯度算子，m为曲率模值，g(·)为依赖于图像的单调递减的扩散函数，I(x,y,t)为t时刻第x行第y列的像素值、构成平滑图像， 为像素值I(x,y,t)对时间t的导数，I(x,y,0)为原始图像，I0(x,y)为0时刻原始图像的像素值；
[0100] 其中，扩散函数表达式采用 x为检测算子，k为阈值。
[0101] 4)用最小均方算法进行阈值寻优，进一步控制扩散函数的扩散强度；
[0102] 根据最小均方算法中的梯度下降法来求取最优阈值k，表达式为式(3)，[0103] k(n+1)＝k(n)-μ▽k(n) (3)
[0104] 式(3)中，k(n+1)为下一时刻的阈值，k(n)为当前时刻的阈值，μ为步长，▽k(n)为均方误差的梯度；
[0105] 其中，根据最小均方算法中的梯度下降法来求取最优阈值k，具体为，[0106] 4-1)对噪声图像进行J层离散小波变换，得到原始信号的小波系数，记为向量Y，且Y＝[AJ,HJ,VJ,DJ,HJ-1,VJ-1,DJ-1,…,H1,V1,D1]，其中，J为层数，向量A为低频系数，向量H为水平高频系数，向量D为对角高频系数；
[0107] 设g(Y)是Y的函数，表达式为式(4)，
[0108]
[0109] 式(4)中， 为基于向量Y对信号w的估计值，故g(Y)为RN到RN的映射，RN为N维向量空间；
[0110] 4-2)基于Sure无偏估计进行阈值的选择，通过风险函数R(t)式(5)来定义，[0111]
[0112] 式(5)中，f、分别为原始信号和原始信号的估计，N为维数；
[0113] 由于小波变换具有正交性，故风险函数可同样在小波域中表达成式(6)，[0114]
[0115] 式(6)中，X、 分别为原始信号的小波系数和原始信号小波系数的估计；
[0116] 4-3)现记Y的风险函数为式(7)，
[0117]
[0118] 则获得其数学期望为式(8)，
[0119]
[0120] 式(8)中，标准差，V为向量Y与原始信号的小波系数X的差值；
[0121] 4-4)可知，V＝Y-X，当V服从高斯分布时，如下式(9)成立，
[0122]
[0123] 式(9)中，Yi为各尺度下的小波系数向量，Vi为各尺度下的Y向量与X向量的差值，i＝1,2,…,N；
[0124] 则得到式(10)，
[0125]
[0126]
[0127] 式(9)中，ti为t在各尺度下的值；
[0128] 4-5)将式(10)代入式(8)，得到风险函数的数学期望ER(t)，为式(11)，[0129]
[0130] 根据Sure无偏估计式(11)，可得均方误差 为式(12)，
[0131]
[0132] 式(12)中，为阈值函数作用于Y后得到的信号估计值，gi为各尺度下 与Y的差值；
[0133] 则均方误差 的梯度表达式为式(13)，
[0134]
[0135] 因为 所以得gi为式(14)，
[0136] gi＝X(kn,Yi)-Yi (14)
[0137] 4-6)将式(14)代入式(13)，得式(15)
[0138]
[0139] 根据式(15)计算出▽k(n)，再根据式(3)计算出k(n+1)，从而获得最优阈值。
[0140] 5)用建立的高保真各向异性滤波模型对提取的高频部分进行处理，对处理后的高频系数和原来的低频系数进行小波重构，得到去噪后图像并输出；
[0141] 将平行于y方向的边缘，忽略y方向的影响和小波变换且Ix＞0，则反映局部信息的检测算子简化为 其中，u＝|W×(▽I+m2)|为检测算子；
[0142] 此时高保真各向异性滤波模型的简化模型为式(16)，
[0143]
[0144] 式(16)中，
[0145]
[0146] 对于式(17)，当在噪声图像的拐点时，由Marr边缘位置知，噪声图像的二阶微分量Ixx＝0，此时有
[0147]
[0148] 又因为 求导得 则得到
[0149]
[0150] 将式(19)代入式(18)中，则得到式(20)，
[0151]
[0152] 对于式(17)，当处于噪声图像的尖峰、角点时，Ix＝0，此时，反映局部信息的检测算子 则式(17)简化为式(21)，
[0153]
[0154] 式(19)、式(20)、式(21)中，k为步骤4)寻优到的最优阈值。
[0155] 6)输出图像，优选为在输出图像之前进行所有像素都直通的判断，若是、则直接输出图像，若否、则返回步骤1)。
[0156] 本发明提供的阈值寻优的高保真各向异性滤波方法，通过用小波变换提取图像的高频部分，在高频部分用二阶微分量曲率模值来反映局部信息，并建立高保真各向异性滤波模型，从而可避免将图像的尖峰、角点误认为噪声，保护了细节信息；并用最小均方算法进行阈值寻优，进一步控制扩散强度，可提高边缘检测的准确性，增强图像滤波效果；用建立的高保真各向异性滤波模型对提取的高频部分进行处理，对处理后的高频系数和原来的低频系数进行小波重构，得到去噪后图像作为输出，去噪性能较现有技术中的经典模型更具优越性，而且运行时间有了明显提高。
[0157] 其中，在高频部分用二阶微分量曲率模值来反映局部信息，曲率模值可反映某点周围曲面变化情况。在图像平滑区域，周围曲面变化细微，曲率模值较小；在图像边缘区域，周围曲面变化较大，曲率模值变得很大。因此，若该点为平坦区域的噪声点时，曲率模值不会如局部方差那样变化很大，所以本发明采用的曲率模值可以较好地区分平坦区域和边缘。
[0158] 现有技术中，经典模型有PM模型(非线性各向异性扩散)和ROF模型(全变分)等，以及最近一些效果较好的模型如NLM模型(非局部平均)，BM3D模型(
块匹配三维协同)，Non-iterative PPB模型(非迭代概率块)，Iterative PPB模型(迭代概率块)，WHT模型(小波硬阈值)，WST模型(小波软阈值)等多种。
[0159] 分析PM模型，通过式(19)，可得其对应简化形式为式(22)，
[0160]
[0161] 比较式(22)与(20)，可得式(20)小于式(22)，即本发明的高保真各向异性滤波模型简化模型小于PM简化模型；所以，当Ix＜k2时， 本发明的高保真各向异性滤波模型使Ix较快地衰减。
[0162] 比较式(22)与式(21)，由于Ixx≥0，则It≥0，故本发明的高保真各向异性滤波模型可抑制灰度值的减小。
[0163] 所以在图像的平坦区域，本发明模型比PM模型能更快地平滑噪声，在图像的边缘纹理等细节信息处，本发明模型虽然有一定的平滑作用，但能更好地保护图像的角点、尖峰等细节信息，本发明模型不仅可有效去除图像噪声，而且保持了图像真实感。
[0164] 仿真实例：
[0165] 用MATLAB
软件对本发明方法进行仿真，得到去噪后图像和数值结果，数值结果可用来评价方法。
[0166] 为验证各滤波算法的有效性和合理性，即采用不同模型进行图像处理，采用像素均为512×512的Nuist图像和Lena图像，如图2所示；分别加方差σ2为20和50的高斯随机噪声，用Matlab软件实行仿真，比较各滤波模型的峰值信噪比(peak signal to noise ratio，PSNR)和结构相似度(structural similarity，SSIM)；其中，SSIM的取值范围为(0,1)，而结构相似度值越接近1表示滤波效果越好。仿真结果得到图像比较为图3-图6。
[0167] 仿真实验对两图像采用不同模型滤波的PSNR与SSIM指标的比较，如下表1所示。
[0168]
[0169] 表1
[0170] 表1中所示的本案模型即可本发明提供的高保真各向异性滤波模型。
[0171] 由表1可知，尽管采用NLM模型处理方法取得了较好的去噪效果，但是对原图像结构信息保护不够，并且从图4d、图6d可以清楚的看到NLM模型处理后的图像丢失了很多边缘纹理等细节信息。
[0172] 如图3-图6所示，从图3c、图5c可以看出，PM模型效果一般图像的一些细节信息较模糊，比如图4c、图6c中，Nuist模型得到Nuist图像的地面、石柱和Lena图像的帽子、眼睛周围等丢失了一些细节信息，这是由于PM模型用了易受噪声影响的梯度来进行边缘检测，且在图像平滑过程中，图像的角点、尖峰与边缘点有相近的梯度模值，若与边缘点以同方式进行平滑，尖峰、角点等细节信息会被磨光。
[0173] 由表1和图3e、图5e可知，BM3D去噪方法不仅有较高的信噪比，视觉效果也较好，但从图4e、图6e可以看出该方法在保护图像细节信息方面有待提高。
[0174] 由表1可知，Iterative PPB方法对图像的处理效果较好，但是该方法耗时过长，不利于实际应用，而Non-iterative PPB方法虽然解决了Iterative PPB方法耗时过长的问题，但从图4f、图6f可以看到，该方法破坏了图像的边缘纹理等细节信息。
[0175] 从图3h、图3i、图5h、图5i、图4h、图4i、图6h、图6i可以看出，WHT方法和WST方法滤波效果不太理想，虽然该方法能够较好地估计噪声方差，并去除图像中的噪声，但会将图像高频子带中的小波系数误认为噪声系数而去除，导致图像的边缘、纹理等细节信息丢失，所以图像存在较严重的“阶梯效应”。
[0176] 然而，从图3j、图5j和图4j、图6j可以看出，采用本发明方法的可视效果最好，因为本发明方法用小波对图像进行分解，对高频部分用二阶微分量进行边缘检测，并用最小均方算法设计阈值进一步控制扩散强度，控制图像的整体结构，能够有效地进行滤波。由表1的评价指标可知，本发明所提方法效果最好，与滤波结果的可视性相一致。
[0177] 由图7可知，本发明提供的高保真各向异性滤波模型的运行时间虽然较经典的PM模型和小波阈值的运行时间慢，但是较NLM方法、BM3D方法、Non-iterative方法和Iterative方法的运行时间快。一方面由于本文用二阶微分量曲率模值进行检测算子，并用最小均方算法进行阈值寻优，能快速地进行各向异性滤波，减少了运行时间；另一方面由于本发明模型采用AOS算法进行数值分解，并用Thomas方法进行求解，从而进一步减少了运行时间。
[0178] 为进一步检测本发明所提方法的性能，在不同噪声方差σ2的条件下，用峰值信噪比对去噪结果进行性能分析，实验结果如图8。由图8可知，本发明所提方法在仿真实验所有方法中具有最高的峰值信噪比，再次证实了本发明所提方法的去噪性能佳。
[0179] 以上显示和描述了本发明的基本原理、主要特征及优点。本行业的技术人员应该了解，本发明不受上述
实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的
权利要求书及其等效物界定。
