function [M] = AllRelativePatchNormXCorr(I,patchsize,range);
 ks = patchsize * patchsize * size(I,3);
 nOri = size(I,3);
 rangerange = 2*range + 1;
 s = size(I);
 D = zeros(size(I,1), size(I,2), nOri);
 n =0;
 I = padimage(I,ceil(patchsize/2),'symmetric');
 LocalMean = zeros(size(I,1),size(I,2),nOri);
 for iOri = 1:nOri
   LocalMean(:,:,iOri) = conv2(ones(patchsize,1),ones(1,patchsize),I(:,:,iOri),'same');
   DENOM_ONESIDE(:,:,iOri) = conv2(ones(patchsize,1),ones(1,patchsize),I(:,:,iOri).*I(:,:,iOri),'same');
 end
 LocalMean = sum(LocalMean,3);
 I = unpadimage(I,ceil(patchsize/2));
 LocalMean = unpadimage(LocalMean,ceil(patchsize/2));
 LocalMean = LocalMean / ks;
 DENOM_ONESIDE = sum(DENOM_ONESIDE,3);
 DENOM_ONESIDE = unpadimage(DENOM_ONESIDE,ceil(patchsize/2));
 DENOM_ONESIDE = DENOM_ONESIDE - ks .* (LocalMean.^2);
 MAX_OF_DENOM = max(DENOM_ONESIDE(:));
 M = zeros(size(I,1),size(I,2),rangerange,rangerange);
 I = padimage(I, patchsize, 'symmetric');
 LocalMean = padimage(LocalMean, patchsize, 'symmetric');
 DENOM_ONESIDE = padimage(DENOM_ONESIDE,patchsize,'symmetric');
 if(any(DENOM_ONESIDE(:) <= 0))
   DENOM_ONESIDE = abs(DENOM_ONESIDE);
   DENOM_ONESIDE = DENOM_ONESIDE + eps;
 end  
 for i = 1:rangerange,
   transi = i - range - 1;
   % for j = floor(rangerange/2):rangerange,
   for j = 1:rangerange,
     n = n+1;
     transj = j - range - 1;
     % if((transj == 0) && (transi == 0)) , continue;, end; % (null translations are too unstable)
     SHIFTI = TranslateIm(I,[transi,transj],.5);
     SHIFTMean = TranslateIm(LocalMean,[transi,transj],.5);
     % NUMERATOR = SHIFTI .* I - LocalMean .* SHIFTI - I .* SHIFTMean + LocalMean .* SHIFTMean;
     NUMERATOR = SHIFTI .* I;
     SHIFTDenom = TranslateIm(DENOM_ONESIDE,[transi,transj],MAX_OF_DENOM);
     for iOri = 1:nOri
       N(:,:,iOri) = conv2(ones(patchsize,1),ones(1,patchsize),NUMERATOR(:,:,iOri),'same');
     end
     N_3 = sum(N,3);
     N_3 = N_3 - ks * sum(LocalMean.*SHIFTMean,3);
     N_3 = RemoveBoundarySignal(N_3,transi,transj,patchsize,0);
     D_3 = (sqrt(DENOM_ONESIDE .* SHIFTDenom));
     D_3 = sum(D_3,3);
     D_3 = RemoveBoundarySignal(D_3,transi,transj,patchsize,MAX_OF_DENOM);
     
     if(0)
      LowerBoundDenom = percentile(D_3(:),.25,1);
      if(LowerBoundDenom == 0)
         LowerBoundDenom = 1/100 * max(D_3(:));
      end
      D_3 = max(D_3,LowerBoundDenom);
     end
     M(:,:,i,j) = unpadimage(N_3 ./ D_3,patchsize);
     % z = M(:,:,i,j);
     % z = z(:);
     % if(any(z<-1))
     %  keyboard
     % end
     % if(i == 13)&&(j ==13), keyboard;,end
   end
 end
% for i = 1:rangerange,
%   transi = i - range - 1;
%   % if(transi == 0), continue;, end; % (null translations are too unstable)
%   for j = 1:(floor(rangerange/2)-1),
%     transj = j - range - 1;
%     % if(transj == 0), continue;, end; % (null translations are too unstable)
%     M(:,:,i,j) = TranslateIm(M(:,:,rangerange-i+1,rangerange-j+1),[transi, transj],0);
%   end
% end


function N_3 = RemoveBoundarySignal(N_3,transi,transj,patchsize,fillval);
sz = size(N_3);
for i = (transi - patchsize):(transi + patchsize)
  for x = 1:sz(2)
    y = mod( (i-1),sz(1))+1;
    N_3(y,x) = fillval;
  end  
end
for j = (transj - patchsize):(transj + patchsize)
  for y = 1:sz(1);
    x = mod( (j-1),sz(2))+1;
    N_3(y,x) = fillval;
  end
end
 
