function  pos_arr   =  Block_matching_RGB(im, Patch_size, Window_size, Group_size, Sliding_step)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Original Author: Weisheng Dong
% Reference: 
%   W. Dong, G. Shi, X. Li, et al., compressive sensing via nonlocal
%   low-rank regularization, IEEE Trans. on Image Processing, vol. 23, no. 8, pp. 3618-3632, 2014.
%
% Modified by Kan Chang, 12/23/2014
% pos_arr = zeros(Group_size, N1*M1 );
% where N1*M1 is the number of the total patches
% note: searching_step = 1 for denosing, 
%       searching_step = 2 for distinguishing sampling and interpolation pixels
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

searching_step = 1;
S         =   Window_size;         % window size, pixel by pixel search
f         =   Patch_size;          % patch size
f2        =   f^2;
s         =   Sliding_step;        % step between two neighbor target patches (to recovery)

N         =   size(im,1)-f+1;      % all the patches in an image (including both the patches to recovery and the patches only for grouping)
M         =   size(im,2)-f+1;
r         =   [1:s:N];             % the patches to recovery

if r(end)+searching_step <= N
    r         =   [r r(end)+searching_step:searching_step:N];     
end
    c         =   [1:s:M];
if c(end)+searching_step <= M
    c         =   [c c(end)+searching_step:searching_step:M];
end
L         =   N*M;
X_R         =   zeros(f*f, L, 'single');
X_G         =   zeros(f*f, L, 'single');
X_B         =   zeros(f*f, L, 'single');

k    =  0;
for i  = 1:f   % collecting simial patches. note that the collecting is pixel by pixel!!!
    for j  = 1:f
        k    =  k+1;
        blk  =  im(i:end-f+i,j:end-f+j,:);
        X_R(k,:) =  reshape(blk(:,:,1),L,1)';
        X_G(k,:) =  reshape(blk(:,:,2),L,1)';
        X_B(k,:) =  reshape(blk(:,:,3),L,1)';
    end
end

% X: patchsize*number of patches(pixel by pixel move)

% Index image
I     =   (1:L);
I     =   reshape(I, N, M);
N1    =   length(r);
M1    =   length(c);
pos_arr   =  zeros(Group_size, N1*M1 );
X_R         =  X_R';
X_G         =  X_G';
X_B         =  X_B';

% searching just go over the selected patches! 
for  i  =  1 : N1       % row
    for  j  =  1 : M1     % column
        
        row     =   r(i);             % i and j are patch position, row and col are pixel position 
        col     =   c(j);
        off     =  (col-1)*N + row;   % column by column scanning, pixel position of current patch
        off1    =  (j-1)*N1 + i;      % column by column scanning, patch position of current patch
                
        rmin    =   max( row-S, 1 );
        rmax    =   min( row+S, N );
        cmin    =   max( col-S, 1 );
        cmax    =   min( col+S, M );
         
        idx     =   I(rmin:searching_step:rmax, cmin:searching_step:cmax);   % searching_step is very important!!!!
        idx     =   idx(:);
        B_R     =   X_R(idx, :);     % B is the pathes that need to be compared within a window
        B_G     =   X_G(idx, :);
        B_B     =   X_B(idx, :);      
        
        v_R     =   X_R(off, :);     % v is the center of the window
        v_G     =   X_G(off, :);
        v_B     =   X_B(off, :);
        
        dis     =   (B_R(:,1) - v_R(1)).^2;
        dis     =   dis + (B_G(:,1) - v_G(1)).^2;
        dis     =   dis + (B_B(:,1) - v_B(1)).^2;
        
        for k = 2:f2
            dis   =  dis + (B_R(:,k) - v_R(k)).^2;
            dis   =  dis + (B_G(:,k) - v_G(k)).^2;
            dis   =  dis + (B_B(:,k) - v_B(k)).^2;
        end
        dis   =  dis./f2;
        [~,ind]   =  sort(dis);        
        pos_arr(:,off1)  =  idx( ind(1:Group_size) );   
        
    end
end