一、小波滤波：
随着小波理论的日益完善，其以自身良好的时频特性在图像去噪领域受到越来越多的关注，开辟了用非线性方法去噪的先河。具体来说，小波能够去噪主要得益于小波变换有如下特点：
（1）低熵性。小波系数的稀疏分布，使图像变换后的熵降低。 意思是对信号（即图像）进行分解后，有更多小波基系数趋于0（噪声），而信号主要部分多集中于某些小波基，采用阈值去噪可以更好的保留原始信号。
（2）多分辨率特性。由于采用了多分辨方法，所以可以非常好地刻画信号的非平稳性，如突变和断点等（例如0-1突变是傅里叶变化无法合理表示的），可以在不同分辨率下根据信号和噪声的分布来消除噪声。
（3）去相关性。小波变换可对信号去相关，且噪声在变换后有白化趋势，所以小波域比时域更利于去噪。
（4）基函数选择灵活。小波变换可灵活选择基函数，也可根据信号特点和去噪要求选择多带小波和小波包等（小波包对高频信号再次分解，可提高时频分辨率），对不同场合，选择不同小波基函数。
根据基于小波系数处理方式的不同，常见去噪方法可分为三类：
（1）基于小波变换模极大值去噪（信号与噪声模极大值在小波变换下会呈现不同变化趋势）
（2）基于相邻尺度小波系数相关性去噪（噪声在小波变换的各尺度间无明显相关性，信号则相反）
（3）基于小波变换阈值去噪
小波去噪实现步骤：
（1）二维信号的小波分解。选择一个小波和小波分解的层次N，然后计算信号s到第N层的分解。
（2）对高频系数进行阈值量化。对于从1~N的每一层，选择一个阈值，并对这一层的高频系数进行软阈值量化处理。
（3）二维小波重构。根据小波分解的第N层的低频系数和经过修改的从第一层到第N的各层高频系数，计算二维信号的小波重构
小波去噪的MATLAB程序：
clc
clear;
zyy = imread('y.jpg'); %读取原图像
subplot(221);
imshow(zyy);
title('原图');
X=imread('chuzao.jpg'); %具有噪声的图片
X=rgb2gray(X);
subplot(222);
imshow(X);
title('待处理噪声');
X=double(X);
%用小波函数coif2对图像X进行2层
% 分解
[c,l]=wavedec2(X,2,'coif2');
% 设置尺度向量
n=[1,2];
% 设置阈值向量 , 对高频小波系数进行阈值处理
p=[10.28,24.08];
nc=wthcoef2('h',c,l,n,p,'s');
% 图像的二维小波重构
X1=waverec2(nc,l,'coif2');
subplot(223);
imshow(uint8(X1));
%colormap(map);
title(' 小波第一次消噪后的图像 ');
%再次对高频小波系数进行阈值处理
mc=wthcoef2('v',nc,l,n,p,'s');
% 图像的二维小波重构
X2=waverec2(mc,l,'coif2');
subplot(224);
imshow(uint8(X2));
title(' 小波重构第二次消噪后的图像 '); 小波去噪的结果：
matlab工作区可以看到运行后的相关数据：
psnr、ssim的值：
在运行完程序之后，在命令窗口输入psnr、ssim口令，如下图所示
%由于X,X1,X2的处理数据为double类型，需要将它们变成uint8类型，才能调用psnr与ssim口令
psnr(zyy,uint8(X))
psnr(zyy,uint8(X1))
psnr(zyy,uint8(X2))
ssim(zyy,uint8(X))
ssim(zyy,uint8(X1))
ssim(zyy,uint8(X2))
%实验结果如下：
%20.7648
%21.5064
%22.3159
%0.3504
%0.3691
%0.3979结果分析：
小波除噪理论上，会得到比较好的除噪效果，但是本次除噪相对均值、中值、高斯等方法除噪效果相对较差；主要原因，是不清楚噪声是高频噪声还是低频噪声，这样会导致部分有用信号丢失；
二、非均值局部滤波(NL-means)：
非局部均值（NL-means）是近年来提出的一项新型的去噪技术。该方法充分利用了图像中的冗余信息，在去噪的同时能最大程度地保持图像的细节特征。基本思想是：当前像素的估计值由图像中与它具有相似邻域结构的像素加权平均得到。
非局部均值去噪算法其实很简单，该种去噪方法和高斯去噪和双边滤波器去噪很像，都是利用一些准则，通过“周围”的像素点加权估计像素点的真实值，如下图所示：
最左边一副图表示Gauss滤波的特点，就是利用图像像素点相近的程度来估计权重，中间幅图表示双边滤波器在考虑像素点本身取值的相近性以外，还考虑了相近像素点与被估计的像素点的距离，如果离被估计的像素点越近将具有更高的权重，非局部均值则是在一个窗口中搜索相近的图像块来进行权重分配（如绿色的框内的区域为搜索的区域，黄色的窗口为搜索的图像块，深褐色的窗口中的中心点为需要去噪的点，通过加权的形式将最相近的几个像素块中的中心点结合起来估计真实值。
MATALB代码：
clc
clear
g = rgb2gray(imread('chuzao.jpg'));
zyy = imread('y.jpg');
g1 = double(g);
%进行NL-means除噪
[m,n]=size(g1);
ds=2;% block size for calculate weight
Ds=5;% search block
h=10;% decay factor
offset=ds+Ds;
PaddedImg = padarray(g1,[ds+Ds,ds+Ds],'symmetric','both');% 扩展图像，便于处理
%距离加权核
%非均值核
[x,y]=meshgrid(-ds:ds,-ds:ds);
kernel=1./(x.*x+y.*y+1);
%均值核
% kernel=ones(2*ds+1,2*ds+1);
kernel=kernel./((2*ds+1)*(2*ds+1));
dst=zeros(m,n);% output
%---------------------------%
% Non-Local Means Denoising
%---------------------------%
for i=1:m
for j=1:n
%当前点坐标和邻域窗口
i1=i+offset;
j1=j+offset;
W1=PaddedImg(i1-ds:i1+ds,j1-ds:j1+ds);
%加权因子矩阵和图像
weight=zeros(2*Ds+1,2*Ds+1);
image=PaddedImg(i1-Ds:i1+Ds,j1-Ds:j1+Ds);
for r=-Ds:Ds
for s=-Ds:Ds
%跳过当前点
if(r==0&&s==0)
continue;
end
%待加权点坐标和邻域窗口
i2=i1+r;
j2=j1+s;
W2=PaddedImg(i2-ds:i2+ds,j2-ds:j2+ds);
%核加权的距离和加权因子
distance=sum(sum(kernel.*(W1-W2).*(W1-W2)));
weight(r+Ds+1,s+Ds+1)=exp(-distance/(h*h));
end
end
%最大权重赋给当前点，归一化权重
weight(Ds+1,Ds+1)=max(max(weight));
weight=weight/(sum(weight(:)));
dst(i,j)=sum(sum(image.*weight));
end
end
%---------------------------%
% output
%---------------------------%
subplot(121),imshow(g1,[]),title('待处理');
subplot(122),imshow(dst,[]),title('NL-means除噪');
psnr(zyy,uint8(g1))
psnr(zyy,uint8(dst))
ssim(zyy,uint8(g1))
ssim(zyy,uint8(dst))
除噪结果：
待除噪的psnr,ssim分别是20.7648 0.3504
除噪后的psnr,ssim分别是26.8054 0.6385
结论分析：
通过结果可以看出NL-means的除噪效果是相当不错的，其它方法psnr与ssim对比可以查看https://blog.csdn.net/A496608119/article/details/110688516
三、BM3D算法除噪
一种非局部去噪方法Non-local method[1]，可以归类到spatial method中，另外用的比较多的还有transform method，基于transform method的方法在image denoise中也取得了很好的效果，不过理论阐述会比较繁琐，如BLS-GSM-Wavelet。
NLM去噪算法使用的是inter-patchcorrelation，而Wavelet shrinkage使用的则是intra-patch correlation。这两种方法都取得了不错的效果，一个很自然的想法就是：可以同时使用他们两个方法吗？这便导出了BM3D去噪算法[2]，算是现在公认的去噪效果最好的算法。
具体内容查看此博主的文章，我没能复现出这个算法，但是通过这个算法，我尝试了将NL-MEAN与小波的融合
四、NL-mean与小波的融合
吸取了BM3D算法的想法，我先对带有噪声的图片进行了小波变换，再对重构的图片进行NL-mean除噪，进而最大化的去除图像噪声；
MATLAB代码：
g = rgb2gray(imread('chuzao.jpg'));
zyy = imread('y.jpg');
%---------------------------%
% 小波除噪
%---------------------------%
g1 = double(g);
%用小波函数coif2对图像X进行2层
% 分解
[c,l]=wavedec2(g1,2,'coif2');
% 设置尺度向量
n=[1,2];
% 设置阈值向量 , 对高频小波系数进行阈值处理
p=[10.28,24.08];
nc=wthcoef2('h',c,l,n,p,'s');
% 图像的二维小波重构
X1=waverec2(nc,l,'coif2');
%再次对高频小波系数进行阈值处理
mc=wthcoef2('v',nc,l,n,p,'s');
% 图像的二维小波重构
g1=waverec2(mc,l,'coif2');
%---------------------------%
% 小波
%---------------------------%
%进行NL-means除噪
[m,n]=size(g1);
ds=2;% block size for calculate weight
Ds=5;% search block
h=10;% decay factor
offset=ds+Ds;
PaddedImg = padarray(g1,[ds+Ds,ds+Ds],'symmetric','both');% 扩展图像，便于处理
%距离加权核
%非均值核
[x,y]=meshgrid(-ds:ds,-ds:ds);
kernel=1./(x.*x+y.*y+1);
%均值核
% kernel=ones(2*ds+1,2*ds+1);
kernel=kernel./((2*ds+1)*(2*ds+1));
dst=zeros(m,n);% output
%---------------------------%
% Non-Local Means Denoising
%---------------------------%
for i=1:m
for j=1:n
%当前点坐标和邻域窗口
i1=i+offset;
j1=j+offset;
W1=PaddedImg(i1-ds:i1+ds,j1-ds:j1+ds);
%加权因子矩阵和图像
weight=zeros(2*Ds+1,2*Ds+1);
image=PaddedImg(i1-Ds:i1+Ds,j1-Ds:j1+Ds);
for r=-Ds:Ds
for s=-Ds:Ds
%跳过当前点
if(r==0&&s==0)
continue;
end
%待加权点坐标和邻域窗口
i2=i1+r;
j2=j1+s;
W2=PaddedImg(i2-ds:i2+ds,j2-ds:j2+ds);
%核加权的距离和加权因子
distance=sum(sum(kernel.*(W1-W2).*(W1-W2)));
weight(r+Ds+1,s+Ds+1)=exp(-distance/(h*h));
end
end
%最大权重赋给当前点，归一化权重
weight(Ds+1,Ds+1)=max(max(weight));
weight=weight/(sum(weight(:)));
dst(i,j)=sum(sum(image.*weight));
end
end
%---------------------------%
% output
%---------------------------%
% subplot(121),imshow(g1,[]),title('均值除噪');
% subplot(122),imshow(dst,[]),title('均值+NL-means除噪');
% subplot(121),imshow(g1,[]),title('中值除噪');
% subplot(122),imshow(dst,[]),title('中值+NL-means除噪');
% subplot(121),imshow(g1,[]),title('高斯除噪');
% subplot(122),imshow(dst,[]),title('高斯+NL-means除噪');
% subplot(121),imshow(g1,[]),title('维纳除噪');
% subplot(122),imshow(dst,[]),title('维纳+NL-means除噪');
subplot(121),imshow(g1,[]),title('小波除噪');
subplot(122),imshow(dst,[]),title('小波+NL-means除噪');
psnr(zyy,uint8(g1))
psnr(zyy,uint8(dst))
ssim(zyy,uint8(g1))
ssim(zyy,uint8(dst))结果：
其它方法混合结果：
psnr与ssim对比：
用不同的滤波方法处理过后的PSNR与SSIM如表一表二所示：（表一的结果，可以通过这篇文章得到https://blog.csdn.net/A496608119/article/details/110688516）
表1 单独滤波方法的PSNR与SSIM
|
|
待除噪图 |
均值滤波 |
中值除噪 |
高斯除噪 |
维纳除噪 |
小波除噪 |
NL-除噪 |
| PSNR |
20.7648 |
23.724 |
24.2671 |
24.0838 |
26.0139 |
22.3159 |
26.8054 |
| SSIM |
0.3504 |
0.5863 |
0.5518 |
0.6006 |
0.7006 |
0.3979 |
0.6385 |
表2 混合滤波方法的 PSNR与SSIM
|
|
待除噪图 |
均值除噪+NL |
中值除噪+NL |
高斯除噪+NL |
维纳除噪+NL |
小波除噪+NL |
| PSNR |
20.7648 |
24.2465 |
25.4840 |
25.838 |
26.4284 |
27.4859 |
| SSIM |
0.3504 |
0.7280 |
0.7695 |
0.7482 |
0.7651 |
0.7308 |
结果分析：
从表一与表二的结果可以看出，混合滤波普遍除噪效果得到提升，尤其是维纳+NL以及小波+NL的效果最明显；这说明了BM3D图像除噪的方法思路是正确的，且除噪效果相对不错；
不足之处：直接将两种方法混合，而不考虑具体内容，虽然psnr与ssim都提高了，但是细节部分变的更模糊了；大家可以根据原理，进一步的优化混合除噪方法；