function [ima_patchs, ima_normalization] = ...
                                   spatial_patchization(im_ini,Patch_width)
% spatial_patchization computes the 3D  collection of patches from the
% noisy images. Patches have size Patch_width x Patch_width, but are 
% coded as vectors.                 
%
%   INPUT:
%     im_ini                 : noisy image
%  	  Patch_width            : patch width 
%
%   OUTPUT:
%     ima_patchs             : 3D collection of patches (each patch is a 
%                              vector)
%    
%     ima_normalization      : "image" coding the number of patches a pixel
%                               pixel belong to.
%
%
%[ima_patchs, ima_normalization] = ...
%                                  spatial_patchization(im_ini,Patch_width)
%
% computes the 3D  collection of patches from the noisy images. Patches 
% have size Patch_width x Patch_width, but are coded as vectors.    
%   
%   Copyright (C) 2012 NL-PCA project
%   Joseph Salmon, Charles-Alban Deledalle, Rebecca Willet, Zachary Harmany
%
%   See The GNU Public License (GPL)

%---------------------------------------------------------------------
%
%   This file is part of NL-PCA.
%
%   NL-PCA is free software: you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as
%   published by the Free Software Foundation, either version 3 of
%   the License, or (at your option) any later version.
%
%   NL-PCA is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public
%   License along with NL-PCA.  If not, see
%   <http://www.gnu.org/licenses/>.
%   Joseph Salmon, Charles-Alban Deledalle, Rebecca Willet, Zachary Harmany
%
%   See The GNU Public License (GPL)
[m,n] = size(im_ini);
ima_patchs=zeros(m-Patch_width+1,n-Patch_width+1,Patch_width^2);
ima_normalization=zeros(m,n);

for i=1:(m-Patch_width+1)
    for j=1:(n-Patch_width+1)
        xrange = mod((i:i+Patch_width-1)-1,m)+1;
        yrange = mod((j:j+Patch_width-1)-1,n)+1;
        B = im_ini(xrange, yrange);
        ima_patchs(i,j,:) = B(:);
        ima_normalization(xrange, yrange) = ...
            ima_normalization(xrange, yrange) + 1;
    end
end
