function [x,y,dx,dy,pkh] = disp_on_blocks_fast(im1,im2,blocksize,overlap,padding,window,method,th,N, Ni);

% DISP_ON_BLOCKS_FAST is the core of the PIV routine that finds the
% displacement field. It computes the displacements of IM2 from IM1.
% The size of the blocks (real resolution of the PIV routine) is determined by the 
% parameter BLOCKSIZE. OVERLAP determines the overlap between different
% blocks. The size of each block of IM1 and IM2 that will be crosscorrelated is 
% (BLOCKSIZE * BLOCKSIZE). IM1 and IM2 can be
% windowed using XWINDOW2D with the method indicated in WINDOW. TH is the
% thershold to determine the centroid of the crosscorrelation peak (out of order). The
% peak is computed on a window of size 2N+1. METHOD is the method used to
% compute the centroid. The process to determine the peak of the
% cross-correlation function is iterative and uses a continous window shift
% (see Gui and Wereley, 2002)
% x and y are the coordinates
% dx is a matrix of x-displacements, dy a matrix of y-displacements, and
% pkh is the height of the croscorrelation peak.
% All units are pixels
% by Xavier Trepat 03/2008
%xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

padding=32; % padding for the edges of the the files /changed to 32 from 16/

sz = size(im1); % size of the image

im1 = xExpandMatrix(im1, 1, 1, padding, padding, padding, padding, 0); % pad
im2 = xExpandMatrix(im2, 1, 1, padding, padding, padding, padding, 0);

inc = round(blocksize*(1-overlap)); % define increment
if (inc<1 || inc>blocksize), error('Wrong Overlap in Correlation Algorithm'), end;

x =   zeros((sz(1)/blocksize-1)*(1/(1-overlap))+1); % position in x
y =   zeros((sz(1)/blocksize-1)*(1/(1-overlap))+1); % position in y
dx =  zeros((sz(1)/blocksize-1)*(1/(1-overlap))+1); % displacements in x
dy =  zeros((sz(1)/blocksize-1)*(1/(1-overlap))+1); % displacements in y
pkh = zeros((sz(1)/blocksize-1)*(1/(1-overlap))+1); % height of the xcorr peak

% MAJOR LOOP
niterations=0;
for ki = 1:inc:sz(1)-blocksize+1,
   fprintf('.');
      for kj =1:inc:sz(1)-blocksize+1,
        niterations=0; DX=inf; DY=inf; % initialize iterative process
        while(niterations<=Ni && abs(dx((ki+inc-1)/inc, (kj+inc-1)/inc ) -DX)>0.02 && abs(dy((ki+inc-1)/inc, (kj+inc-1)/inc ) -DY)>0.02); % set iteration conditions
           
            niterations=niterations+1;
            DX=(dx((ki+inc-1)/inc, (kj+inc-1)/inc ));
            DY=(dy((ki+inc-1)/inc, (kj+inc-1)/inc ));
            
            im11 = im1(ki : ki+blocksize+2*padding-1 , kj : kj+blocksize+2*padding-1); % crop the block with padding
            im22 = im2(ki : ki+blocksize+2*padding-1 , kj : kj+blocksize+2*padding-1);
            
            im11=im11/mean2(im11);
            im22=im22/mean2(im22);
            
            if(DX || DY) %skip first iteration
            im11 = xsubpix_shift(im11,DX/2,DY/2); % subpixel shift of the image
            im22 = xsubpix_shift(im22,-DX/2,-DY/2);
            end
            
            im11 = im11(padding+1 : padding+blocksize , padding+1 : padding+blocksize); % crop the block
            im22 = im22(padding+1 : padding+blocksize , padding+1 : padding+blocksize);
            
            im11 = im11 - mean2(im11); % subtract the mean
            im22 = im22 - mean2(im22);                  
            
            if(window),   
                im11 = xwindow2D(im11,window); % multiply each block by a window
                im22 = xwindow2D(im22,window);                
            end
            
            c = xcorrf2(im22,im11) / (std(im11(:))*std(im22(:))*blocksize^2); % compute the correlation function in Fourier Space and normalize
            c = real(c(1:end-1,1:end-1)); % resize for x_cntr 
            
            x((ki+inc-1)/inc, (kj+inc-1)/inc ) = kj+blocksize/2;
            y((ki+inc-1)/inc, (kj+inc-1)/inc ) = ki+blocksize/2;         
            dx((ki+inc-1)/inc, (kj+inc-1)/inc ) = DX + x_cntr(c,N,method,'X',th); % compute the displacement
            dy((ki+inc-1)/inc, (kj+inc-1)/inc ) = DY + x_cntr(c,N,method,'Y',th);   
            pkh((ki+inc-1)/inc, (kj+inc-1)/inc ) = max2(c); % store peak height
            
        end
    end
end
fprintf('\n');
