传统的高斯滤波,均值滤波,为局部滤波,即对周围邻域的点加权生成当前点,加权因子反应出周围点对当前点的影响,这些加权因子基于某种理论获得,如高斯滤波基于低通,均值滤波认为点与点之间的影响是均匀的。

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)
<https://blog.csdn.net/u010839382/article/details/48241929>