non-local Means(非局部均值)降噪算法及快速算法原理與
Non-Local Means算法原理：
Non-Local Means顧名思義，這是一種非局部平均算法。何為局部平均濾波算法呢？那是在一個目標像素周圍區域平滑取均值的方法，所以非局部均值濾波就意味着它使用圖像中的所有像素，這些像素根據某種相似度進行加權平均。濾波后圖像清晰度高，而且不丟失細節。
非局部均值濾波由Baudes提出，其出發點應該是借鑒了越多幅圖像加權的效果越好的現象，那么在同一幅圖像中對具有相同性質的區域進行分類並加權平均得到去噪后的圖片，應該降噪效果也會越好。該算法使用自然圖像中普遍存在的冗余信息來去噪聲。與雙線性濾波、中值濾波等利用圖像局部信息來濾波不同，它利用了整幅圖像進行去噪。即以圖像塊為單位在圖像中尋找相似區域，再對這些區域取平均，較好地濾除圖像中的高斯噪聲。NL-Means的濾波過程可以用下面公式來表示：
w代表權重。衡量相似度的方法有很多，最常用的是根據兩個像素亮度差值的平方來估計。由於有噪聲，單獨的一個像素並不可靠，所以使用它們的鄰域，只有鄰域相似度高才能說這兩個像素的相似度高。衡量兩個圖像塊的相似度最常用的方法是計算他們之間的歐氏距離：
其中 a 是高斯核的標准差。在求歐式距離的時候，不同位置的像素的權重是不一樣的，距離塊的中心越近，權重越大，距離中心越遠，權重越小，權重服從高斯分布。實際計算中考慮到計算量的問題，常常采用均勻分布的權重。
如上圖所示，p為去噪的點，q1和q2的鄰域與p相似，所以權重w(p,q1) 和w(p,q2) 較大，而鄰域相差比較大的點q3的權重值w(p,q3) 很小。如果用一幅圖把所有點的權重表示出來，那就得到下面這些權重圖：
一般而言，考慮到算法復雜度，搜索區域大概取21x21，相似度比較的塊的可以取7x7。實際中，常常需要根據噪聲來選取合適的參數。當高斯噪聲的標准差σ 越大時，為了提高算法魯棒性，需要增大塊區域，同樣也需要增加搜索區域。同時，濾波系數h 與 σ 滿足正相關：h=kσ，當塊變大時，k 需要適當減小。
Non-Local Means算法實現：
function [output]=NLmeansfilter(input,t,f,h)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% input: image to be filtered
% t: radio of search window
% f: radio of similarity window
% h: degree of filtering
%
% Author: Jose Vicente Manjon Herrera & Antoni Buades
% Date: 09-03-2006
%
% Implementation of the Non local filter proposed for A. Buades, B. Coll and J.M. Morel in
% "A non-local algorithm for image denoising"
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Size of the image
[m n]=size(input);
% Memory for the output
Output=zeros(m,n);
% Replicate the boundaries of the input image
input2 = padarray(input,[f f],'symmetric');
% Used kernel
kernel = make_kernel(f);
kernel = kernel / sum(sum(kernel));
h=h*h;
for i=1:m
for j=1:n
i1 = i+ f;
j1 = j+ f;
W1= input2(i1-f:i1+f , j1-f:j1+f);
wmax=0;
average=0;
sweight=0;
rmin = max(i1-t,f+1);
rmax = min(i1+t,m+f);
smin = max(j1-t,f+1);
smax = min(j1+t,n+f);
for r=rmin:1:rmax
for s=smin:1:smax
if(r==i1 && s==j1) continue; end;
W2= input2(r-f:r+f , s-f:s+f);
d = sum(sum(kernel.*(W1-W2).*(W1-W2)));
w=exp(-d/h);
if w>wmax
wmax=w;
end
sweight = sweight + w;
average = average + w*input2(r,s);
end
end
average = average + wmax*input2(i1,j1);
sweight = sweight + wmax;
if sweight > 0
output(i,j) = average / sweight;
else
output(i,j) = input(i,j);
end
end
end
function [kernel] = make_kernel(f)
kernel=zeros(2*f+1,2*f+1);
for d=1:f
value= 1 / (2*d+1)^2 ;
for i=-d:d
for j=-d:d
kernel(f+1-i,f+1-j)= kernel(f+1-i,f+1-j) + value ;
end
end
end
kernel = kernel ./ f;
下面是我寫的測試函數：
%% 測試函數
clc,clear all,close all;
ima=double(imread('3.jpg'));
[wid,len,channels]=size(ima);
% add noise
sigma=10;
rima=ima+sigma*randn(size(ima));
% denoise
fima=rima;
if channels>2
fori=1:channels
fima(:,:,i)=NLmeansfilter(rima(:,:,i),5,2,sigma);
end
end
% show results
subplot(1,3,1),imshow(uint8(ima)),title('original');
subplot(1,3,2),imshow(uint8(rima)),title('noisy');
subplot(1,3,3),imshow(uint8(fima)),title('filtered');
由於原始算法的復雜度較高，導致算法耗時及較長，所以針對NLM算法產生了不少優化算法，如使用積分圖像技術對算法進行加速。為了降低空間復雜度，將偏移量作為最外層循環，即每次只需要在一個偏移方向上求取積分圖像，並對該積分圖像進行處理。而不需要一次性求取出所有積分圖像，參考【6】。算法流程見下圖：
先構造一個關於像素差值的積分圖像：
其中
這樣在計算兩個鄰域和 間的距離時，就可以在常量時間內完成：
這樣，整個算法復雜度將降為 。具體代碼如下：
function DenoisedImg=fastNLmeans(I,ds,Ds,h)
%I:含噪聲圖像
%ds:鄰域窗口半徑
%Ds:搜索窗口半徑
%h:高斯函數平滑參數
%DenoisedImg：去噪圖像
I=double(I);
[m,n]=size(I);
PaddedImg = padarray(I,[Ds+ds+1,Ds+ds+1],'symmetric','both');
PaddedV = padarray(I,[Ds,Ds],'symmetric','both');
average=zeros(m,n);
sweight=average;
wmax=average;
h2=h*h;
d2=(2*ds+1)^2;
for t1=-Ds:Ds
for t2=-Ds:Ds
if(t1==0&&t2==0)
continue;
end
St=integralImgSqDiff(PaddedImg,Ds,t1,t2);
v = PaddedV(1+Ds+t1:end-Ds+t1,1+Ds+t2:end-Ds+t2);
w=zeros(m,n);
for i=1:m
for j=1:n
i1=i+ds+1;
j1=j+ds+1;
Dist2=St(i1+ds,j1+ds)+St(i1-ds-1,j1-ds-1)-St(i1+ds,j1-ds-1)-St(i1-ds-1,j1+ds);
Dist2=Dist2/d2;
w(i,j)=exp(-Dist2/h2);
sweight(i,j)=sweight(i,j)+w(i,j);
average(i,j)=average(i,j)+w(i,j)*v(i,j);
end
end
wmax=max(wmax,w);
end
end
average=average+wmax.*I;
sweight=sweight+wmax;
DenoisedImg=average./sweight;
function Sd = integralImgSqDiff(PaddedImg,Ds,t1,t2)
%PaddedImg：邊緣填充后的圖像
%Ds:搜索窗口半徑
%(t1,t2):偏移量
%Sd:積分圖像
[m,n]=size(PaddedImg);
m1=m-2*Ds;
n1=n-2*Ds;
Sd=zeros(m1,n1);
Dist2=(PaddedImg(1+Ds:end-Ds,1+Ds:end-Ds)-PaddedImg(1+Ds+t1:end-Ds+t1,1+Ds+t2:end-Ds+t2)).^2;
for i=1:m1
for j=1:n1
if i==1 && j==1
Sd(i,j)=Dist2(i,j);
elseif i==1 && j~=1
Sd(i,j)=Sd(i,j-1)+Dist2(i,j);
elseif i~=1 && j==1
Sd(i,j)=Sd(i-1,j)+Dist2(i,j);
else
Sd(i,j)=Dist2(i,j)+Sd(i-1,j)+Sd(i,j-1)-Sd(i-1,j-1);
end
end
end
方案2：
function DenoisedImg=fastNLmeans2(I,ds,Ds,h)
I=double(I);
[m,n]=size(I);
PaddedImg = padarray(I,[Ds+ds+1,Ds+ds+1],'symmetric','both');
PaddedV = padarray(I,[Ds,Ds],'symmetric','both');
average=zeros(m,n);
wmax=average;
sweight=average;
h2=h*h;
d=(2*ds+1)^2;
for t1=-Ds:Ds
for t2=-Ds:Ds
if(t1==0&&t2==0)
continue;
end
Sd=integralImgSqDiff(PaddedImg,Ds,t1,t2);
SqDist2=Sd(2*ds+2:end-1,2*ds+2:end-1)+Sd(1:end-2*ds-2,1:end-2*ds-2)...
-Sd(2*ds+2:end-1,1:end-2*ds-2)-Sd(1:end-2*ds-2,2*ds+2:end-1);
SqDist2=SqDist2/d;
w=exp(-SqDist2/h2);
v = PaddedV(1+Ds+t1:end-Ds+t1,1+Ds+t2:end-Ds+t2);
average=average+w.*v;
wmax=max(wmax,w);
sweight=sweight+w;
end
end
average=average+wmax.*I;
average=average./(wmax+sweight);
DenoisedImg = average;
function Sd = integralImgSqDiff(PaddedImg,Ds,t1,t2)
Dist2=(PaddedImg(1+Ds:end-Ds,1+Ds:end-Ds)-PaddedImg(1+Ds+t1:end-Ds+t1,1+Ds+t2:end-Ds+t2)).^2;
Sd = cumsum(Dist2,1);
Sd = cumsum(Sd,2);
測試結果如下：
參考：
《a non-local algorithm for image denoising》[J].IEEE
https://en.wikipedia.org/wiki/Non-local_means
http://wenhuix.github.io/research/denoise.html
http://blog.csdn.net/piaoxuezhong/article/details/78317861
http://blog.csdn.net/tuyang120428941/article/details/7052487
http://blog.csdn.net/u010839382/article/details/48241929
————————————————
版權聲明：本文為CSDN博主「Naruto_Q」的原創文章，遵循CC 4.0 BY-SA版權協議，轉載請附上原文出處鏈接及本聲明。
原文鏈接：https://blog.csdn.net/piaoxuezhong/article/details/78345929