function ima = reprojection_UWA(ima_patchs,xg,yg)

%   reprojection_UWA compute the reprojection step that from a collections 
%   of patches representing an image, provide an image estimate. The 
%   various pixels estimates (due to overlapping) are uniformly average.
%                 
%   INPUT:
%     ima_patchs             : denoised collections of patches 
%  	  xg,yg                  : potential subpart of the image
%
%   OUTPUT:
%    ima                     : denoised images by uniformly averaging the
%                              patches containing a pixel of interest to
%                              proposed a denoised pixel estimator
%
%   ima = reprojection_UWA(ima_patchs,xg,yg)
%   
%   Compute the reprojection step that from a collections 
%   of patches ima_patchs representing an image, provide an image estimate.
%   The various pixels estimates (due to overlapping) are uniformly
%   average.  Potentially restrict the process to a smaller window,
%   centered around xg,yg
%   
%   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, P] = size(ima_patchs);
    if nargin <= 2
        xg = zeros(M,N);
        yg = zeros(M,N);
    end
    p = sqrt(P);
    ima = zeros(M+p-1, N+p-1);
    norm = zeros(M+p-1, N+p-1);
    for i = 1:M
        for j = 1:N
            xrange = mod(round((i-1)+(1:p)+xg(i,j))-1,M+p-1)+1;
            yrange = mod(round((j-1)+(1:p)+yg(i,j))-1,N+p-1)+1;
            ima(xrange, yrange) = ...
                ima(xrange, yrange) + ...
                reshape(ima_patchs(i,j,:), p, p);
            norm(xrange, yrange) = ...
                norm(xrange, yrange) + ...
                1;
        end
    end
    ima = ima ./ norm;
    clear norm;
