CN105913382B - 阈值寻优的高保真各向异性滤波方法 - Google Patents
阈值寻优的高保真各向异性滤波方法 Download PDFInfo
- Publication number
- CN105913382B CN105913382B CN201610115932.XA CN201610115932A CN105913382B CN 105913382 B CN105913382 B CN 105913382B CN 201610115932 A CN201610115932 A CN 201610115932A CN 105913382 B CN105913382 B CN 105913382B
- Authority
- CN
- China
- Prior art keywords
- formula
- image
- fidelity
- wavelet
- equation
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000001914 filtration Methods 0.000 title claims abstract description 55
- 238000000034 method Methods 0.000 title claims abstract description 53
- 238000009792 diffusion process Methods 0.000 claims abstract description 25
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 20
- 230000000694 effects Effects 0.000 claims abstract description 18
- 230000009466 transformation Effects 0.000 claims abstract description 11
- 239000013598 vector Substances 0.000 claims description 29
- 238000001514 detection method Methods 0.000 claims description 16
- 238000012545 processing Methods 0.000 claims description 14
- 150000001875 compounds Chemical class 0.000 claims description 9
- 238000005457 optimization Methods 0.000 claims description 8
- 238000011478 gradient descent method Methods 0.000 claims description 6
- 239000011541 reaction mixture Substances 0.000 claims description 6
- 238000000354 decomposition reaction Methods 0.000 claims description 5
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 claims description 4
- 230000003247 decreasing effect Effects 0.000 claims description 3
- 238000007781 pre-processing Methods 0.000 claims description 3
- 238000006467 substitution reaction Methods 0.000 claims description 3
- 230000001419 dependent effect Effects 0.000 claims 1
- 238000003708 edge detection Methods 0.000 abstract description 9
- 230000009286 beneficial effect Effects 0.000 abstract description 4
- 238000000605 extraction Methods 0.000 abstract description 3
- 238000004088 simulation Methods 0.000 description 6
- 230000008569 process Effects 0.000 description 3
- 238000004458 analytical method Methods 0.000 description 2
- 230000008859 change Effects 0.000 description 2
- 230000002093 peripheral effect Effects 0.000 description 2
- 230000009467 reduction Effects 0.000 description 2
- 230000000007 visual effect Effects 0.000 description 2
- 238000007476 Maximum Likelihood Methods 0.000 description 1
- 238000004364 calculation method Methods 0.000 description 1
- 238000004891 communication Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000011156 evaluation Methods 0.000 description 1
- 238000009499 grossing Methods 0.000 description 1
- 238000003706 image smoothing Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000003672 processing method Methods 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 230000011218 segmentation Effects 0.000 description 1
- 239000004575 stone Substances 0.000 description 1
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
- G06T5/70—Denoising; Smoothing
Landscapes
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- Theoretical Computer Science (AREA)
- Image Processing (AREA)
- Image Analysis (AREA)
Abstract
本发明公开了一种阈值寻优的高保真各向异性滤波方法，首先用小波变换提取图像的高频部分，在高频部分用二阶微分量曲率模值来反映局部信息，并建立高保真各向异性滤波模型；再用最小均方算法进行阈值寻优，进一步控制扩散强度；然后用建立的高保真各向异性滤波模型对提取的高频部分进行处理，对处理后的高频系数和原来的低频系数进行小波重构，得到去噪后图像作为输出。有益效果：可避免将图像的尖峰、角点误认为噪声，保护了细节信息，还可提高边缘检测的准确性，增强图像滤波效果，故去噪性能较现有技术中的经典模型更具优越性，而且运行时间有了明显提高。
Description
技术领域
本发明涉及一种图像去噪方法，特别是涉及一种阈值寻优的高保真各向异性滤波方法，属于图像处理技术领域。
背景技术
图像去噪是图像处理与计算机领域中最基本的问题之一。图像去噪是图像的边缘检测、增强和分割等图像处理中的一个非常重要预处理。现阶段，图像去噪方法种类繁多，有基于偏微分方程、概率论、小波等图像去噪方法。
近20年来，基于偏微分方程的方法在图像去噪中得到了广泛的应用。Perona和Malik首次提出了经典的各向异性扩散模型，即PM模型，该模型由一个关于图像梯度模值的扩散函数控制扩散强度；自PM模型以后，基于PDE的去噪方法继续发展，其有效地提高了边缘检测能力。
目前，Zhang K B，Gao X B等人提出了一种NLM滤波器，该方法较好地保持了图像纹理和细节特征，但是计算复杂高，处理速度慢。BM3D去噪方法主要通过块匹配以及三维变换域滤波等技术进行两次去噪，第一次基础去噪为第二次最终去噪提供权值，使得最终去噪效果更加好，该方法去噪后的图像不仅有较高的信噪比，视觉效果也很好。Deledalle等提出了PPB滤波器，该滤波器通过重新定义像素之间的相似度，并用加权最大似然估计完成图像的重建，在合成孔径雷达(SAR)图像的处理上达到较好的效果，但实验结果表明该滤波器运行时间较长，而非迭代滤波方法不能保留图像纹理等细节特征，SAR图像噪声的先验信息未充分利用。
小波具有良好的时频变换特性，使得它在图像去噪领域发挥着重要作用，在基于小波变换的去噪方法中，小波硬阈值去噪算法和小波软阈值去噪算法是较为有效地降噪手段。目前，基于小波的图像阈值去噪取得了很多研究成果，在扩散滤波以及它的相应变化方法方面，人们逐渐对离散小波变换和偏微分方程之间的联系产生了兴趣。
在各向异性扩散滤波方法中，边缘检测的准确性对滤波结果有重要的影响。然而，目前基于偏微分方程去噪方法的扩散强度大都使用梯度信息检测边缘，在图像的尖峰、窄边缘处一阶偏导数为0，对图像边缘纹理的尖峰进行了平滑处理，磨平了图像的尖峰、角点等细节信息，从而使图像失真。
总之，上述传统方法时效性低，复杂度较高，且在处理过程中图像降噪不稳定，有明显的“阶梯”效应，并且保边缘性效果不是很好。
发明内容
本发明的主要目的在于，克服现有技术中的不足，提供一种阈值寻优的高保真各向异性滤波方法，可有效保护图像的尖峰、角点等细节信息，提高边缘检测的准确性，增强图像滤波效果，去噪性能优异；而且运行时间短，有利于实际应用。
为了达到上述目的，本发明所采用的技术方案是：
一种阈值寻优的高保真各向异性滤波方法，包括以下步骤：
1)输入图像，采用各向异性扩散模型、即PM模型对图像进行预处理后输出含噪声的噪声图像；
2)判断是否考虑局部特征处理，若是、则进入步骤3)，若否、则进入步骤6)；
3)用小波变换提取噪声图像的高频部分，在高频部分用二阶微分量曲率模值m来反映局部信息，并建立高保真各向异性滤波模型；
根据式(1)计算曲率模值m，
式(1)中，Ixx、Ixy分别为Ix在x,y方向的方向导数，Iyy为Iy在y方向的方向导数，其中，Ix为图像I在x方向分量，Iy为图像I在y方向分量；
高保真各向异性滤波模型为式(2)，
式(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时刻原始图像的像素值；
其中，扩散函数表达式采用x为检测算子，k为阈值；
4)用最小均方算法进行阈值寻优，进一步控制扩散函数的扩散强度；
根据最小均方算法中的梯度下降法来求取最优阈值k，表达式为式(3)，
k(n+1)＝k(n)-μ▽k(n) (3)
式(3)中，k(n+1)为下一时刻的阈值，k(n)为当前时刻的阈值，μ为步长，▽k(n)为均方误差的梯度；
5)用建立的高保真各向异性滤波模型对提取的高频部分进行处理，对处理后的高频系数和原来的低频系数进行小波重构，得到去噪后图像并输出；
将平行于y方向的边缘，忽略y方向的影响和小波变换且Ix＞0，则反映局部信息的检测算子简化为其中，u＝|W×(▽I+m2)|为检测算子；
此时高保真各向异性滤波模型的简化模型为式(16)，
式(16)中，
6)输出图像。
本发明进一步设置为：所述步骤4)根据最小均方算法中的梯度下降法来求取最优阈值k，具体为，
4-1)对噪声图像进行J层离散小波变换，得到原始信号的小波系数，记为向量Y，且Y＝[AJ,HJ,VJ,DJ,HJ-1,VJ-1,DJ-1,…,H1,V1,D1]，其中，J为层数，向量A为低频系数，向量H为水平高频系数，向量D为对角高频系数；
设g(Y)是Y的函数，表达式为式(4)，
式(4)中，为基于向量Y对信号w的估计值，故g(Y)为RN到RN的映射，RN为N维向量空间；
4-2)基于Sure无偏估计进行阈值的选择，通过风险函数R(t)式(5)来定义，
式(5)中，f、分别为原始信号和原始信号的估计，N为维数；
由于小波变换具有正交性，故风险函数可同样在小波域中表达成式(6)，
式(6)中，X、分别为原始信号的小波系数和原始信号小波系数的估计；
4-3)现记Y的风险函数为式(7)，
则获得其数学期望为式(8)，
式(8)中，σn为标准差，V为向量Y与原始信号的小波系数X的差值；
4-4)可知，V＝Y-X，当V服从高斯分布时，如下式(9)成立，
式(9)中，Yi为各尺度下的小波系数向量，Vi为各尺度下的Y向量与X向量的差值，i＝1,2,…,N；
则得到式(10)，
式(9)中，ti为t在各尺度下的值；
4-5)将式(10)代入式(8)，得到风险函数的数学期望ER(t)，为式(11)，
根据Sure无偏估计式(11)，可得均方误差为式(12)，
式(12)中，为阈值函数作用于Y后得到的信号估计值，gi为各尺度下与Y的差值；
则均方误差的梯度表达式为式(13)，
因为所以得gi为式(14)，
gi＝X(kn,Yi)-Yi (14)
4-6)将式(14)代入式(13)，得式(15)
根据式(15)计算出▽k(n)，再根据式(3)计算出k(n+1)，从而获得最优阈值。
本发明进一步设置为：所述步骤5)中的式(17)，当在噪声图像的拐点时，由Marr边缘位置知，噪声图像的二阶微分量Ixx＝0，此时有
又因为求导得则得到
将式(19)代入式(18)中，则得到式(20)，
式(19)、式(20)中，k为步骤4)寻优到的最优阈值。
本发明进一步设置为：所述步骤5)中的式(17)，当处于噪声图像的尖峰、角点时，Ix＝0，此时，反映局部信息的检测算子则式(17)简化为式(21)，
式(21)中，k为步骤4)寻优到的最优阈值。
本发明进一步设置为：所述步骤6)输出图像之前进行所有像素都直通的判断，若是、则直接输出图像，若否、则返回步骤1)。
与现有技术相比，本发明具有的有益效果是：
1、在去噪性能方面，利用二阶微分量曲率模值反映图像局部信息，可避免将图像的尖峰、角点等误认为是噪声，有效地保护了图像的尖峰、角点等细节信息，用最小均方算法设计阈值，提高了边缘检测的准确性，增强了图像滤波的效果，去噪性能更具优越性。
2、在运行时间方面，通过用二阶微分量曲率模值进行检测算子，并用最小均方算法进行阈值寻优，能快速地进行各向异性滤波，从而减少运行时间，而且又采用AOS算法进行数值分解，并用Thomas方法进行求解，可进一步减少运行时间。
上述内容仅是本发明技术方案的概述，为了更清楚的了解本发明的技术手段，下面结合附图对本发明作进一步的描述。
附图说明
图1为本发明阈值寻优的高保真各向异性滤波方法的流程图；
图2为本发明仿真实验所用的原始图像，依次为Nuist图像和Lena图像；
图3为原始图像分别采用不同模型的滤波结果(方差σ2＝20)；
图4为原始图像分别采用不同模型滤波的边缘提取(方差σ2＝20)；
图5为原始图像分别采用不同模型的滤波结果(方差σ2＝50)；
图6为原始图像分别采用不同模型滤波的边缘提取(方差σ2＝50)；
图7为原始图像分别采用不同模型滤波的运行时间；
图8为原始图像分别采用不同模型在不同方差噪声下峰值信噪比。
具体实施方式
下面结合说明书附图，对本发明作进一步的说明。
如图1所示，本发明提供一种阈值寻优的高保真各向异性滤波方法，包括以下步骤：
1)输入图像，采用各向异性扩散模型、即PM模型对图像进行预处理后输出含噪声的噪声图像。
2)判断是否考虑局部特征处理，若是、则进入步骤3)，若否、则进入步骤6)。
3)用小波变换提取噪声图像的高频部分，在高频部分用二阶微分量曲率模值m来反映局部信息，并建立高保真各向异性滤波模型；
根据式(1)计算曲率模值m，
式(1)中，Ixx、Ixy分别为Ix在x,y方向的方向导数，Iyy为Iy在y方向的方向导数，其中，Ix为图像I在x方向分量，Iy为图像I在y方向分量；
高保真各向异性滤波模型为式(2)，
式(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时刻原始图像的像素值；
其中，扩散函数表达式采用x为检测算子，k为阈值。
4)用最小均方算法进行阈值寻优，进一步控制扩散函数的扩散强度；
根据最小均方算法中的梯度下降法来求取最优阈值k，表达式为式(3)，
k(n+1)＝k(n)-μ▽k(n) (3)
式(3)中，k(n+1)为下一时刻的阈值，k(n)为当前时刻的阈值，μ为步长，▽k(n)为均方误差的梯度；
其中，根据最小均方算法中的梯度下降法来求取最优阈值k，具体为，
4-1)对噪声图像进行J层离散小波变换，得到原始信号的小波系数，记为向量Y，且Y＝[AJ,HJ,VJ,DJ,HJ-1,VJ-1,DJ-1,…,H1,V1,D1]，其中，J为层数，向量A为低频系数，向量H为水平高频系数，向量D为对角高频系数；
设g(Y)是Y的函数，表达式为式(4)，
式(4)中，为基于向量Y对信号w的估计值，故g(Y)为RN到RN的映射，RN为N维向量空间；
4-2)基于Sure无偏估计进行阈值的选择，通过风险函数R(t)式(5)来定义，
式(5)中，f、分别为原始信号和原始信号的估计，N为维数；
由于小波变换具有正交性，故风险函数可同样在小波域中表达成式(6)，
式(6)中，X、分别为原始信号的小波系数和原始信号小波系数的估计；
4-3)现记Y的风险函数为式(7)，
则获得其数学期望为式(8)，
式(8)中，标准差，V为向量Y与原始信号的小波系数X的差值；
4-4)可知，V＝Y-X，当V服从高斯分布时，如下式(9)成立，
式(9)中，Yi为各尺度下的小波系数向量，Vi为各尺度下的Y向量与X向量的差值，i＝1,2,…,N；
则得到式(10)，
式(9)中，ti为t在各尺度下的值；
4-5)将式(10)代入式(8)，得到风险函数的数学期望ER(t)，为式(11)，
根据Sure无偏估计式(11)，可得均方误差为式(12)，
式(12)中，为阈值函数作用于Y后得到的信号估计值，gi为各尺度下与Y的差值；
则均方误差的梯度表达式为式(13)，
因为所以得gi为式(14)，
gi＝X(kn,Yi)-Yi (14)
4-6)将式(14)代入式(13)，得式(15)
根据式(15)计算出▽k(n)，再根据式(3)计算出k(n+1)，从而获得最优阈值。
5)用建立的高保真各向异性滤波模型对提取的高频部分进行处理，对处理后的高频系数和原来的低频系数进行小波重构，得到去噪后图像并输出；
将平行于y方向的边缘，忽略y方向的影响和小波变换且Ix＞0，则反映局部信息的检测算子简化为其中，u＝|W×(▽I+m2)|为检测算子；
此时高保真各向异性滤波模型的简化模型为式(16)，
式(16)中，
对于式(17)，当在噪声图像的拐点时，由Marr边缘位置知，噪声图像的二阶微分量Ixx＝0，此时有
又因为求导得则得到
将式(19)代入式(18)中，则得到式(20)，
对于式(17)，当处于噪声图像的尖峰、角点时，Ix＝0，此时，反映局部信息的检测算子则式(17)简化为式(21)，
式(19)、式(20)、式(21)中，k为步骤4)寻优到的最优阈值。
6)输出图像，优选为在输出图像之前进行所有像素都直通的判断，若是、则直接输出图像，若否、则返回步骤1)。
本发明提供的阈值寻优的高保真各向异性滤波方法，通过用小波变换提取图像的高频部分，在高频部分用二阶微分量曲率模值来反映局部信息，并建立高保真各向异性滤波模型，从而可避免将图像的尖峰、角点误认为噪声，保护了细节信息；并用最小均方算法进行阈值寻优，进一步控制扩散强度，可提高边缘检测的准确性，增强图像滤波效果；用建立的高保真各向异性滤波模型对提取的高频部分进行处理，对处理后的高频系数和原来的低频系数进行小波重构，得到去噪后图像作为输出，去噪性能较现有技术中的经典模型更具优越性，而且运行时间有了明显提高。
其中，在高频部分用二阶微分量曲率模值来反映局部信息，曲率模值可反映某点周围曲面变化情况。在图像平滑区域，周围曲面变化细微，曲率模值较小；在图像边缘区域，周围曲面变化较大，曲率模值变得很大。因此，若该点为平坦区域的噪声点时，曲率模值不会如局部方差那样变化很大，所以本发明采用的曲率模值可以较好地区分平坦区域和边缘。
现有技术中，经典模型有PM模型(非线性各向异性扩散)和ROF模型(全变分)等，以及最近一些效果较好的模型如NLM模型(非局部平均)，BM3D模型(块匹配三维协同)，Non-iterative PPB模型(非迭代概率块)，Iterative PPB模型(迭代概率块)，WHT模型(小波硬阈值)，WST模型(小波软阈值)等多种。
分析PM模型，通过式(19)，可得其对应简化形式为式(22)，
比较式(22)与(20)，可得式(20)小于式(22)，即本发明的高保真各向异性滤波模型简化模型小于PM简化模型；所以，当Ix＜k2时，本发明的高保真各向异性滤波模型使Ix较快地衰减。
比较式(22)与式(21)，由于Ixx≥0，则It≥0，故本发明的高保真各向异性滤波模型可抑制灰度值的减小。
所以在图像的平坦区域，本发明模型比PM模型能更快地平滑噪声，在图像的边缘纹理等细节信息处，本发明模型虽然有一定的平滑作用，但能更好地保护图像的角点、尖峰等细节信息，本发明模型不仅可有效去除图像噪声，而且保持了图像真实感。
仿真实例：
用MATLAB软件对本发明方法进行仿真，得到去噪后图像和数值结果，数值结果可用来评价方法。
为验证各滤波算法的有效性和合理性，即采用不同模型进行图像处理，采用像素均为512×512的Nuist图像和Lena图像，如图2所示；分别加方差σ2为20和50的高斯随机噪声，用Matlab软件实行仿真，比较各滤波模型的峰值信噪比(peak signal to noiseratio，PSNR)和结构相似度(structural similarity，SSIM)；其中，SSIM的取值范围为(0,1)，而结构相似度值越接近1表示滤波效果越好。仿真结果得到图像比较为图3-图6。
仿真实验对两图像采用不同模型滤波的PSNR与SSIM指标的比较，如下表1所示。
表1
表1中所示的本案模型即可本发明提供的高保真各向异性滤波模型。
由表1可知，尽管采用NLM模型处理方法取得了较好的去噪效果，但是对原图像结构信息保护不够，并且从图4d、图6d可以清楚的看到NLM模型处理后的图像丢失了很多边缘纹理等细节信息。
如图3-图6所示，从图3c、图5c可以看出，PM模型效果一般图像的一些细节信息较模糊，比如图4c、图6c中，Nuist模型得到Nuist图像的地面、石柱和Lena图像的帽子、眼睛周围等丢失了一些细节信息，这是由于PM模型用了易受噪声影响的梯度来进行边缘检测，且在图像平滑过程中，图像的角点、尖峰与边缘点有相近的梯度模值，若与边缘点以同方式进行平滑，尖峰、角点等细节信息会被磨光。
由表1和图3e、图5e可知，BM3D去噪方法不仅有较高的信噪比，视觉效果也较好，但从图4e、图6e可以看出该方法在保护图像细节信息方面有待提高。
由表1可知，Iterative PPB方法对图像的处理效果较好，但是该方法耗时过长，不利于实际应用，而Non-iterative PPB方法虽然解决了Iterative PPB方法耗时过长的问题，但从图4f、图6f可以看到，该方法破坏了图像的边缘纹理等细节信息。
从图3h、图3i、图5h、图5i、图4h、图4i、图6h、图6i可以看出，WHT方法和WST方法滤波效果不太理想，虽然该方法能够较好地估计噪声方差，并去除图像中的噪声，但会将图像高频子带中的小波系数误认为噪声系数而去除，导致图像的边缘、纹理等细节信息丢失，所以图像存在较严重的“阶梯效应”。
然而，从图3j、图5j和图4j、图6j可以看出，采用本发明方法的可视效果最好，因为本发明方法用小波对图像进行分解，对高频部分用二阶微分量进行边缘检测，并用最小均方算法设计阈值进一步控制扩散强度，控制图像的整体结构，能够有效地进行滤波。由表1的评价指标可知，本发明所提方法效果最好，与滤波结果的可视性相一致。
由图7可知，本发明提供的高保真各向异性滤波模型的运行时间虽然较经典的PM模型和小波阈值的运行时间慢，但是较NLM方法、BM3D方法、Non-iterative方法和Iterative方法的运行时间快。一方面由于本文用二阶微分量曲率模值进行检测算子，并用最小均方算法进行阈值寻优，能快速地进行各向异性滤波，减少了运行时间；另一方面由于本发明模型采用AOS算法进行数值分解，并用Thomas方法进行求解，从而进一步减少了运行时间。
为进一步检测本发明所提方法的性能，在不同噪声方差σ2的条件下，用峰值信噪比对去噪结果进行性能分析，实验结果如图8。由图8可知，本发明所提方法在仿真实验所有方法中具有最高的峰值信噪比，再次证实了本发明所提方法的去噪性能佳。
以上显示和描述了本发明的基本原理、主要特征及优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书及其等效物界定。
Claims (5)
1.一种阈值寻优的高保真各向异性滤波方法，其特征在于，包括以下步骤：
1)输入图像，采用各向异性扩散模型、即PM模型对图像进行预处理后输出含噪声的噪声图像；
2)判断是否考虑局部特征处理，若是、则进入步骤3)，若否、则进入步骤6)；
3)用小波变换提取噪声图像的高频部分，在高频部分用二阶微分量曲率模值m来反映局部信息，并建立高保真各向异性滤波模型；
根据式(1)计算曲率模值m，
式(1)中，Ixx、Ixy分别为Ix在x,y方向的方向导数，Iyy为Iy在y方向的方向导数，其中，Ix为图像I在x方向分量，Iy为图像I在y方向分量；
高保真各向异性滤波模型为式(2)，
式(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时刻原始图像的像素值；
其中，扩散函数表达式采用x为检测算子，k为阈值；
4)用最小均方算法进行阈值寻优，进一步控制扩散函数的扩散强度；
根据最小均方算法中的梯度下降法来求取最优阈值k，表达式为式(3)，
式(3)中，k(n+1)为下一时刻的阈值，k(n)为当前时刻的阈值，μ为步长，为均方误差的梯度；
5)用建立的高保真各向异性滤波模型对提取的高频部分进行处理，对处理后的高频系数和原来的低频系数进行小波重构，得到去噪后图像并输出；
将平行于y方向的边缘，忽略y方向的影响和小波变换，且Ix＞0，则反映局部信息的检测算子简化为其中，为检测算子；
此时高保真各向异性滤波模型的简化模型为式(16)，
式(16)中，
6)输出图像。
2.根据权利要求1所述的阈值寻优的高保真各向异性滤波方法，其特征在于：所述步骤4)根据最小均方算法中的梯度下降法来求取最优阈值k，具体为，
4-1)对噪声图像进行J层离散小波变换，得到原始信号的小波系数，记为向量Y，且
Y＝[AJ,HJ,ZJ,DJ,HJ-1,ZJ-1,DJ-1,…,H1,Z1,D1]，其中，J为层数，向量A为低频系数，向量H为水平高频系数，向量D为对角高频系数；
设g(Y)是Y的函数，表达式为式(4)，
式(4)中，为基于向量Y对信号w的估计值，故g(Y)为RN到RN的映射，RN为N维向量空间；
4-2)基于Sure无偏估计进行阈值的选择，通过风险函数R(t)式(5)来定义，
式(5)中，f、分别为原始信号和原始信号的估计，N为维数；
由于小波变换具有正交性，故风险函数同样在小波域中表达成式(6)，
式(6)中，X、分别为原始信号的小波系数值和原始信号小波系数的估计；
4-3)现记Y的风险函数为式(7)，
则获得其数学期望为式(8)，
式(8)中，σn为标准差，V为向量Y与原始信号的小波系数值X的差值；
4-4)可知，V＝Y-X，当V服从高斯分布时，如下式(9)成立，
式(9)中，Yi为各尺度下的小波系数向量，Vi为各尺度下的Y向量与X的差值，P为(|Yi|＞t)的概率，i＝1,2,…,N；
则得到式(10)，
式(10)中，ti为t在各尺度下的值；
4-5)将式(10)代入式(8)，得到风险函数的数学期望ER(t)，为式(11)，
根据Sure无偏估计式(11)，可得均方误差为式(12)，
式(12)中，为阈值函数作用于Y后得到的信号估计值，gi为各尺度下与Y的差值；
则均方误差的梯度表达式为式(13)，
因为所以得gi为式(14)，
gi＝X(kn,Yi)-Yi (14)
4-6)将式(14)代入式(13)，得式(15)
根据式(15)计算出再根据式(3)计算出k(n+1)，从而获得最优阈值。
3.根据权利要求1所述的阈值寻优的高保真各向异性滤波方法，其特征在于：所述步骤5)中的式(17)，当在噪声图像的拐点时，由Marr边缘位置知，噪声图像的二阶微分量Ixx＝0，此时有
又因为求导得则得到
将式(19)代入式(18)中，则得到式(20)，
式(19)、式(20)中，k为步骤4)寻优到的最优阈值。
4.根据权利要求1所述的阈值寻优的高保真各向异性滤波方法，其特征在于：所述步骤5)中的式(17)，当处于噪声图像的尖峰、角点时，Ix＝0，此时，反映局部信息的检测算子则式(17)简化为式(21)，
式(21)中，k为步骤4)寻优到的最优阈值。
5.根据权利要求1所述的阈值寻优的高保真各向异性滤波方法，其特征在于：所述步骤6)输出图像之前进行所有像素都直通的判断，若是、则直接输出图像，若否、则返回步骤1)。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201610115932.XA CN105913382B (zh) | 2016-03-01 | 2016-03-01 | 阈值寻优的高保真各向异性滤波方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201610115932.XA CN105913382B (zh) | 2016-03-01 | 2016-03-01 | 阈值寻优的高保真各向异性滤波方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN105913382A CN105913382A (zh) | 2016-08-31 |
| CN105913382B true CN105913382B (zh) | 2018-11-23 |
Family
ID=56744474
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201610115932.XA Active CN105913382B (zh) | 2016-03-01 | 2016-03-01 | 阈值寻优的高保真各向异性滤波方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN105913382B (zh) |
Families Citing this family (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107247819B (zh) * | 2017-05-02 | 2020-07-24 | 歌尔科技有限公司 | 传感器的滤波方法和滤波器 |
| CN108120698B (zh) * | 2017-11-22 | 2020-05-19 | 南京航空航天大学 | 面向柔性薄板结构载荷分布监测的光纤层析成像方法 |
| CN109345465B (zh) * | 2018-08-08 | 2023-04-07 | 西安电子科技大学 | 基于gpu的高分辨率图像实时增强方法 |
| CN109697704A (zh) * | 2018-11-28 | 2019-04-30 | 山东师范大学 | 基于bm3d算法的自适应全变分espi图像降噪方法及系统 |
| CN110009582A (zh) * | 2019-03-28 | 2019-07-12 | 华南理工大学 | 一种基于曲率特征的各向异性图像去噪方法 |
| CN113567603B (zh) * | 2021-07-22 | 2022-09-30 | 华谱科仪(大连)科技有限公司 | 色谱谱图的检测分析方法及电子设备 |
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104463810A (zh) * | 2014-12-25 | 2015-03-25 | 南京信息工程大学 | 基于tv流的自适应扩散滤波图像去噪算法 |
| CN104463811A (zh) * | 2014-12-29 | 2015-03-25 | 南京信息工程大学 | 基于能量泛函的图像平滑与锐化算法 |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104077746B (zh) * | 2013-03-29 | 2017-03-01 | 富士通株式会社 | 灰度图像处理方法及其装置 |
-
2016
- 2016-03-01 CN CN201610115932.XA patent/CN105913382B/zh active Active
Patent Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104463810A (zh) * | 2014-12-25 | 2015-03-25 | 南京信息工程大学 | 基于tv流的自适应扩散滤波图像去噪算法 |
| CN104463811A (zh) * | 2014-12-29 | 2015-03-25 | 南京信息工程大学 | 基于能量泛函的图像平滑与锐化算法 |
Non-Patent Citations (4)
| Title |
|---|
| Robust adaptive directional lifting wavelet transform for image denoising;X.T Wang等;《IET Image Processing》;20111231;第5卷(第3期);第249-260页 * |
| Single Image Super-Resolution With Non-Local Means and Steering Kernel Regression;Kaibing Zhan等;《IEEE Transactions on Image Processing》;20121130;第21卷(第11期);第4544-4556页 * |
| 基于脉冲藕合神经网络和图像嫡的各向异性扩散模型研究;郭业才等;《物理学报》;20151231;第64卷(第19期);第194204-1-194204-11页 * |
| 非线性扩散图像去噪中的藕合自适应保真项研究;朱立新等;《计算机辅助设计与图形学学报》;20061031;第18卷(第10期);第1519-1524页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN105913382A (zh) | 2016-08-31 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN105913382B (zh) | 阈值寻优的高保真各向异性滤波方法 | |
| Zhao et al. | Detail-preserving image denoising via adaptive clustering and progressive PCA thresholding | |
| CN105005975B (zh) | 基于pcnn和图像熵的各向异性扩散的图像去噪方法 | |
| CN103473755B (zh) | 基于变化检测的sar图像稀疏去噪方法 | |
| CN111242854B (zh) | 一种图像去噪方法 | |
| CN111223063A (zh) | 基于纹理特征和双核函数的手指静脉图像nlm去噪方法 | |
| CN108805845A (zh) | 小波包与偏微分方程相结合的图像去噪方法 | |
| CN117726537A (zh) | 自适应多尺度特征融合AMFFD-Net的SAR图像去噪网络方法及系统 | |
| Wu et al. | Pyramid edge detection based on stack filter | |
| CN102314675B (zh) | 基于小波高频的贝叶斯去噪方法 | |
| CN111091107A (zh) | 一种人脸区域边缘检测方法、装置及存储介质 | |
| CN104616266B (zh) | 一种基于广义自回归异方差模型的噪声方差估计方法 | |
| Yu et al. | EPLL image denoising with multi-feature dictionaries | |
| Li et al. | Improved non-local means algorithm for image denosing | |
| Potocnik et al. | Image enhancement by using directional wavelet transform | |
| CN112927169B (zh) | 一种基于小波变换和改进的加权核范数最小化的遥感影像去噪方法 | |
| Gantenapalli et al. | Selective mean filtering for reducing impulse noise in digital color images | |
| CN109840896B (zh) | 一种基于梯度和自适应曲率特征的图像去噪方法 | |
| CN110009591B (zh) | 一种基于四阶偏微分方程的自适应阈值图像去噪方法 | |
| CN113484913A (zh) | 一种多粒度特征融合卷积神经网络的地震资料去噪方法 | |
| CN106651804B (zh) | 基于自适应扩散滤波的车载图像去噪方法及系统 | |
| CN113191961B (zh) | 一种基于双边滤波和特征识别的超声去噪方法及系统 | |
| Han et al. | ICA-domain filtering of Poisson noise images | |
| Li et al. | Understanding neural-network denoisers through an activation function perspective | |
| Shimodaira | Patch-based bilateral filter: Algorithms and performance |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | Publication | ||
| PB01 | Publication | ||
| C10 | Entry into substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant | ||
| CP02 | Change in the address of a patent holder |
Address after: Room 420, block C, Kechuang headquarters building, No. 320, pubin Road, Jiangpu street, Pukou District, Nanjing, Jiangsu Patentee after: Nanjing University of Information Science and Technology Address before: The Olympic Avenue in Jianye District of Nanjing city of Jiangsu Province, No. 69 210019 Patentee before: Nanjing University of Information Science and Technology |
|
| CP02 | Change in the address of a patent holder | ||
| EE01 | Entry into force of recordation of patent licensing contract | ||
| EE01 | Entry into force of recordation of patent licensing contract |
Application publication date: 20160831 Assignee: Wuxi Yunqi Technology Information Co.,Ltd. Assignor: Nanjing University of Information Science and Technology Contract record no.: X2023320000234 Denomination of invention: High fidelity anisotropic filtering method for threshold optimization Granted publication date: 20181123 License type: Common License Record date: 20231121 |
|
| EE01 | Entry into force of recordation of patent licensing contract | ||
| EE01 | Entry into force of recordation of patent licensing contract |
Application publication date: 20160831 Assignee: Nanjing Xinqihang Software Technology Co.,Ltd. Assignor: Nanjing University of Information Science and Technology Contract record no.: X2023980051736 Denomination of invention: High fidelity anisotropic filtering method for threshold optimization Granted publication date: 20181123 License type: Common License Record date: 20231213 Application publication date: 20160831 Assignee: Fujian Kailan Information Technology Co.,Ltd. Assignor: Nanjing University of Information Science and Technology Contract record no.: X2023980051725 Denomination of invention: High fidelity anisotropic filtering method for threshold optimization Granted publication date: 20181123 License type: Common License Record date: 20231213 |
|
| EC01 | Cancellation of recordation of patent licensing contract | ||
| EC01 | Cancellation of recordation of patent licensing contract |
Assignee: Nanjing Xinqihang Software Technology Co.,Ltd. Assignor: Nanjing University of Information Science and Technology Contract record no.: X2023980051736 Date of cancellation: 20240507 Assignee: Fujian Kailan Information Technology Co.,Ltd. Assignor: Nanjing University of Information Science and Technology Contract record no.: X2023980051725 Date of cancellation: 20240507 |