1. 程式人生 > >non-local Means(非區域性均值)降噪演算法及快速演算法原理與實現

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