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