快速塊比對的非局部均值去噪算法_Fast Block Matching Non local means

function  [posIdx, weiIdx]   =  BM_NL_means_test(Img,par)
% Img ： 輸入圖像
% posIdx ： 與中心塊相似的塊的位置索引
% weiIdx ： 與中心塊相似的塊的權重索引

search_r  =  par.s_r; %search block的半徑
similarWinNum =  par.nblk; % 相似window的個數
win_r           =  par.win_r; %similarity window的半徑
win=2*win_r+1;%
winSize     =  win^2; % 每個similarity window的元素個數
hp            =  par.hp^2; % 高斯參數
step=1;
%在圖像四周鏡面填充厚度為win_r的像素點，行數和列數增加了win（2*win_r）個
%填充的目的是為了讓每一個像素都能被計算到
PaddedImg = padarray(Img,[win_r,win_r],'symmetric','both');
[R , C] = size(PaddedImg);

%  R=row+2*win_r C=col+2*win_r
N          =  R - win + 1; % 最後一個window的起始行坐标（步長1）
M          =  C - win + 1; % 最後一個window的起始列坐标
winNum   =  N * M; % window 的總個數

leftUpRow  =  1:step:N; % 每個window的行起始坐标
leftUpRow  =  [leftUpRow leftUpRow(end)+1:N]; % 添加最後一個塊，因為可能最後不一定正好取完
leftUpCol  =  1:step:M; % 每個window的列起始坐标
leftUpCol  =  [leftUpCol leftUpCol(end)+1:M]; % 添加最後一個塊，因為可能最後不一定正好取完


%分塊：按列取一維塊中不同位置的資料（連續的位址，速度可比正常提高10倍）!!!!!!
%這一小段是W dong學者BlockMatching中最精華部分！！！！
X    =  zeros(winSize, M*N, 'single');
k    =  0;
for i  = 1:win
    for j  = 1:win
        k    =  k+1;%每次取R*C個window的第k個像素值
        blk  =  PaddedImg(i:end-win+i,j:end-win+j);%這是算法精華所在噢~
        X(k,:) =  blk(:)';%每列皆是一個window的win*win個像素，共R2*C2列，完成降維
    end
end
X = X';

I      = reshape(1:winNum, N, M); % 所有塊的索引
N1     = length(leftUpRow); % 多少行中心塊
M1     = length(leftUpCol); % 多少列中心塊中心塊
posIdx = zeros(similarWinNum, N1*M1 ); % 每個中心塊的相似塊索引
weiIdx = zeros(similarWinNum, N1*M1 ); % 每個中心塊對應的相似塊的權重
%
rmin   = bsxfun(@max, bsxfun(@minus, leftUpRow, search_r), 1); % 搜尋窗的行最小坐标
rmax   = bsxfun(@min, bsxfun(@plus, leftUpRow, search_r), N); % 搜尋窗的行最大坐标
cmin   = bsxfun(@max, bsxfun(@minus, leftUpCol, search_r), 1); % 搜尋窗的列最小坐标
cmax   = bsxfun(@min, bsxfun(@plus, leftUpCol, search_r), M); % 搜尋窗的列最大坐标
%其實 經過padarry填充以及後來減去win寬度後，此處的N1,N2恰好是初始噪聲圖的次元一模一樣的噢~（R,C）
% step=1是M1/N1是和M/N一樣大的
for  i  =  1 : N1
    for  j  =  1 : M1 
        
        offInAllWin   = (leftUpCol(j)-1)*N + leftUpRow(i); % 目前中心塊在所有塊中的索引
        offInCentralWin  = (j-1)*N1 + i; % 目前中心塊在所有中心塊中的索引

        idx   = I(rmin(i):rmax(i), cmin(j):cmax(j));  %在由[rmin:rmax, cmin:cmax]确定的搜尋窗中搜尋相似塊
        idx   = idx(:);%目前搜尋框中所有塊的索引

        B       =   X(idx, :);        
        v       =   X(offInAllWin, :);  
        %将搜尋視窗的所有window都與目前中心window求它們間的距離
        dis     =   (B(:,1) - v(1)).^2;
        for k = 2:winSize
            dis   =  dis + (B(:,k) - v(k)).^2;
        end
        dis   = dis ./ (winSize); % 歸一化
       
        similarWinInd = minN(dis, similarWinNum); % 找到距離最小的幾個塊的索引
        posIdx(:,offInCentralWin)  =  idx(similarWinInd); % 儲存相似塊索引       
        
        wei   =  exp( -dis(similarWinNum) ./ hp ); % 高斯
        weiIdx(:,offInCentralWin)  =  wei ./ (sum(wei(:))+eps); % 歸一化             
    end
end

function index = minN(data, N)
% 功能 ：找到 data 中最小的 N 個數，并傳回索引
[~, index] = sort(data);
index = index(1:N);
           

 主函數

par.s_r       =5; % search window的半徑
par.nblk      =35; % 用來權重和去噪的相似塊的個數（當然了，也可以是search win所有塊）
par.win_r     =2; % similarity window的半徑
par.step=1;
par.hp       =  13.5; % 高斯平滑參數
tic
I=double(imread('taylor.jpg'));
randn('state',0);%保證每次的高斯随機白噪聲是一樣的随機序列
Img=I+10*randn(size(I));
[row , col] = size(Img);
[posIdx, weiIdx]=BM_NL_means_test(Img,par);
pixel_posIdx = double(Img(posIdx));
Recon_X=sum((pixel_posIdx.*weiIdx),1);
Denoising_Im=reshape(Recon_X,row,col);
toc
           

% figure,imshow(I/256);%原圖
% figure,imshow(Img/256);%加噪聲後
% figure,imshow(Denoising_Im/256);%去噪後
           

參考文獻和資料：

1.http://see.xidian.edu.cn/faculty/wsdong/wsdong_Publication.htm 《Sparsity-based Image Denoising via Dictionary Learning and Structural》源碼中find_block_fast.m

2.《A non-local algorithm for image denoising》Antoni Buades, Bartomeu Coll等人

3. http://blog.csdn.net/wenxuegeng/article/details/51476980 Non Local Means-塊比對MATLAB和GPU實作

4. http://www.cnblogs.com/luo-peng/p/4785922.html 非局部均值去噪（NL-means）

若有疑問，歡迎郵件交流~