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）