Non-Local Means 非局部均值去噪滤
传统的高斯滤波，均值滤波，为局部滤波，即对周围邻域的点加权生成当前点，加权因子反应出周围点对当前点的影响，这些加权因子基于某种理论获得，如高斯滤波基于低通，均值滤波认为点与点之间的影响是均匀的。
1.经典的Non-Local Means 滤波
Non-local Means 非局部均值去噪滤波可以视为局部均值滤波的特例，它的目的是使用与当前点纹理类似的区域，对当前点加权。也即加权因子，是基于被加权点与当前点的邻域的相似性产生，即：
其中I是一个较大范围的搜索/加权框，w(x,y)是依赖邻域【黑灰色部分】算出的权重：
w(x,y)一般定义为一个与欧式距离(2范数)相关的函数，设x，y的邻域宏块的欧式距离为d，即
d=||block(x)-block(y)||/block_size
则y加权到x点的加权因子为
w(x,y)=exp(-(d*d/(h*h)))
h为衰减因子，h越小，加权因子越小，则加权点对当前点的影响越小，一般边缘保持得好但是噪声会严重，反之则边缘保持差图像更加光滑。
实际操作中，要更新当前点，先计算出以当前点为中心的搜索框I所有点的加权因子，取最大的加权因子付给当前点位置，然后对于这个同搜索框尺寸加权矩阵W进行归一化，最终当前点的结果为：
p=sum(W.*I)
计算欧式距离时，有时会考虑周围点对中心点的影响，会利用核函数对欧式距离加权，即加权因子重写为为：
w(x,y)=exp(-(k*d*d/(h*h)))
使用核函数对距离进行加权的matlab代码为：
clc;
clear all;
close all;
%---------------------------%
% input
%---------------------------%
src=imread('test.jpg');
src=rgb2gray(src);
src=double(src);
figure,imshow(src,[]),title('src image')
%---------------------------%
% parameter
%---------------------------%
[m,n]=size(src);
ds=2;% block size for calculate weight
Ds=5;% search block
h=10;% decay factor
offset=ds+Ds;
PaddedImg = padarray(src,[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
%---------------------------%
figure,imshow(dst,[]),title('dst');
效果如下:
2.使用积分图加速
假设图像共像个素点，搜索窗口大小，领域窗口大小， 计算两个矩形邻域间相似度的时间为，对于每个像素点需要计算它与搜索窗口内个像素间的相似度，故NL-means复杂度为
积分图加速本质是，将图像的所有点，计算某一偏离坐标点方向的权重，一次性计算出来，而不是单次计算某一点的所有偏离点。具体，计算当前图像与一定偏移后的图像的diff并平方，继而计算与原图像同等尺寸的积分图，则每一个点与偏离它在一定偏移的坐标点的距离，通过积分图就可以计算出来了。积分图的作用体现在使用线性计算，减少重复计算，即存储换时间。
具体源码如下:
clc;
clear all;
close all;
%---------------------------%
% input
%---------------------------%
src=imread('lena.jpg');
src=rgb2gray(src);
src=double(src);
figure,imshow(src,[]),title('src image')
%---------------------------%
% parameter
%---------------------------%
[m,n]=size(src);
ds=2;% block size for calculate weight
Ds=5;% search block
h=10;% decay factor
offset=ds+Ds;
PaddedImg = padarray(src,[ds+Ds,ds+Ds],'symmetric','both');% 扩展图像，便于处理
%---------------------------%
% Non-Local Means Denoising
%---------------------------%
sumimage=zeros(m,n);
sumweight=zeros(m,n);
maxweight=zeros(m,n);
image=PaddedImg(1+Ds:Ds+m+ds,1+Ds:Ds+n+ds);
[M,N]=size(image);
for r=-Ds:Ds
for s=-Ds:Ds
%跳过当前点偏移
if(r==0&&s==0)
continue;
end
%求得差值积分图
wimage=PaddedImg(1+Ds+r:Ds+m+ds+r,1+Ds+s:Ds+n+ds+s);
diff=image-wimage;
diff=diff.^2;
J=cumsum(diff,1);
J=cumsum(J,2);
%计算距离
distance=J(M-m+1:M,N-n+1:N)+J(1:m,1:n)-J(M-m+1:M,1:n)-J(1:m,N-n+1:N);
distance=distance/((2*ds+1).^2);
%计算权重并获得单个偏移下的加权图像
weight=exp(-distance./(h*h));
sumimage=sumimage+weight.*wimage(ds+1:ds+m,ds+1:ds+n);
sumweight=sumweight+weight;
maxweight=max(maxweight,weight);
end
end
sumimage=sumimage+maxweight.*image(ds+1:ds+m,ds+1:ds+n);
sumweight=sumweight+maxweight;
dst=sumimage./sumweight;
%---------------------------%
% output
%---------------------------%
figure,imshow(dst,[]),title('dst');
效果如下：
参考：
[1]FromentJ. Parameter-Free Fast Pixelwise Non-Local Means Denoising[J]. Image ProcessingOn Line, 2014, 4: 300-326
[2] 非局部均值去噪（NL-means）
————————————————
版权声明：本文为CSDN博主「攻城狮凌风」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
原文链接：https://blog.csdn.net/qianhen123/article/details/81043217