function [x2a,bestshifth,bestshiftv] = imalign_(x2,x1,maxshift,maxshiftchange,delta)
%
% function [x2a,bestshifth,bestshiftv] = imalign_(x2,x1,maxshift,maxshiftchange,delta)
%
% Written by: Luís Borges de Almeida
% email : luis.almeida@lx.it.pt
%
% Aligns image x2 relative to x1.
% maxshift - maximal aligment size (in pixels).
%
% Optional parameters:
%   - maxshiftchange: maxshiftchange = 0 indicates full search of all shift values (from -maxshift to maxshift). 
%                     Maxshiftchange default value: 0.
%   - delta: Shifts are computed for squares of size delta x delta, based on correlations computed on
%            squares of size (2 x delta + 1) x (2 x delta + 1). Delta default value: 100.

[h,v] = size(x1);

if ~exist('delta')
    delta = 100;   % Shifts are computed for squares of size delta x delta, based on correlations computed on
                   % squares of size (2 x delta + 1) x (2 x delta + 1)
end

if ~exist('maxshiftchange')
   maxshiftchange = 0;         % By default we perform a full search of all shift values
end

% myx1 = double(x1);
% myx2 = double(x2);

% Set the means of both images to zero (for internal processing only). We don't need to normalize variances

%myx1 = myx1 - mean(mean(myx1));
%myx2 = myx2 - mean(mean(myx2));

indh = 1;
for ch = 1+delta+maxshift:delta:h-delta-maxshift
    
    fprintf('%d',ch)
    
    indv = 1;
    for cv = 1+delta+maxshift:delta:v-delta-maxshift

        fprintf('.')
        pause(0.01)

        win1 = x1(ch-delta:ch+delta,cv-delta:cv+delta);
        win1 = win1 - mean(mean(win1));
        
        maxval = -1;

        if (indh == 1) | (maxshiftchange == 0)   % Perform full search for shift values, in first row
           
           for shifth = -maxshift:maxshift
              for shiftv = -maxshift:maxshift
                 win2 = x2(ch+shifth-delta:ch+shifth+delta,cv+shiftv-delta:cv+shiftv+delta);
                 win2 = win2 - mean(mean(win2));
                 corr = sum(sum(win1 .* win2)) / sqrt(sum(sum(win2 .* win2)));
                 %corr = sum(sum(win1 .* win2));
                 if corr > maxval
                    maxval = corr;
                    bestshifth(indh,indv) = shifth;
                    bestshiftv(indh,indv) = shiftv;
                 end
              end
           end

        else     % In all other rows, search for previous shift +-maxshiftchange

           for shifth = max(bestshifth(indh-1,indv)-maxshiftchange,-maxshift):min(bestshifth(indh-1,indv)+maxshiftchange,maxshift)
              for shiftv = max(bestshiftv(indh-1,indv)-maxshiftchange,-maxshift):min(bestshiftv(indh-1,indv)+maxshiftchange,maxshift)
                 win2 = x2(ch+shifth-delta:ch+shifth+delta,cv+shiftv-delta:cv+shiftv+delta);
                 win2 = win2 - mean(mean(win2));
                 corr = sum(sum(win1 .* win2)) / sqrt(sum(sum(win2 .* win2)));
                 %corr = sum(sum(win1 .* win2));
                 if corr > maxval
                    maxval = corr;
                    bestshifth(indh,indv) = shifth;
                    bestshiftv(indh,indv) = shiftv;
                 end
              end
           end

        end

           indv = indv + 1;
       end
       indh = indh + 1;
    
    fprintf('\n')
end


% Perform the actual alignment

%x2a = zeros(size(x2));     % This forces x2a to be of type double!

delta2 = ceil(delta / 2);

indh = 1;
for ch = 1+delta+maxshift:delta:h-delta-maxshift
    
    %fprintf('%d',ch)
    
    indv = 1;

    for cv = 1+delta+maxshift:delta:v-delta-maxshift

        %fprintf('.')
        
        shifth = bestshifth(indh,indv);
        shiftv = bestshiftv(indh,indv);

        x2a(ch-delta2:ch+delta2,cv-delta2:cv+delta2) = x2(ch+shifth-delta2:ch+shifth+delta2,cv+shiftv-delta2:cv+shiftv+delta2);
                
        indv = indv + 1;
    end
    
    % Write the end and the beginning of the "column".
    
    x2a(ch-delta2:ch+delta2,cv+delta2+1:v-max(0,shiftv)) = x2(ch+shifth-delta2:ch+shifth+delta2,cv+shiftv+delta2+1:v+min(shiftv,0));
    if shiftv > 0
        for i=shiftv-1:-1:0
            x2a(ch-delta2:ch+delta2,v-i) = x2(ch+shifth-delta2:ch+shifth+delta2,v);
        end
    end
    
    cv1 = 1+delta+maxshift;
    shifth = bestshifth(indh,1);
    shiftv = bestshiftv(indh,1);
    x2a(ch-delta2:ch+delta2,1-min(0,shiftv):cv1-delta2-1) = x2(ch+shifth-delta2:ch+shifth+delta2,1+max(shiftv,0):cv1-delta2-1+shiftv);
    if shiftv<0
        for i=0:-shiftv-1
            x2a(ch-delta2:ch+delta2,1+i) = x2(ch+shifth-delta2:ch+shifth+delta2,1);
        end
    end
    
    indh = indh + 1;
    
    %fprintf('\n')
end

% Write the end and the beginning of the "rows". This can lead to index overflow problems...

indv = 1;
indh = indh - 1;
ch1 = 1+delta+maxshift;
for cv = 1+delta+maxshift:delta:v-delta-maxshift
    shifth = bestshifth(indh,indv);
    shiftv = bestshiftv(indh,indv);
    x2a(ch+delta2+1:h-max(0,shifth),cv-delta2:cv+delta2) = x2(ch+shifth+delta2+1:h+min(shifth,0),cv+shiftv-delta2:cv+shiftv+delta2);
    if shifth>0
        for i=0:shifth-1
            x2a(h-i,cv-delta2:cv+delta2) = x2(h,cv+shiftv-delta2:cv+shiftv+delta2);
        end
    end
            
    
    shifth = bestshifth(1,indv);
    shiftv = bestshiftv(1,indv);
    x2a(1-min(0,shifth):ch1-delta2-1,cv-delta2:cv+delta2) = x2(1+max(shifth,0):ch1-delta2-1+shifth,cv+shiftv-delta2:cv+shiftv+delta2);
    if shifth<0
        for i=0:-shifth-1
            x2a(1+i,cv-delta2:cv+delta2) = x2(1,cv+shiftv-delta2:cv+shiftv+delta2);
        end
    end
            
    
    indv = indv + 1;
    
end

indv = indv - 1;


shifth = bestshifth(1,1);
shiftv = bestshiftv(1,1);
x2a(1-min(0,shifth):ch1-delta2-1,1-min(0,shiftv):cv1-delta2-1) = x2(1+max(shifth,0):ch1-delta2-1+shifth,1+max(shiftv,0):cv1-delta2-1+shiftv);
if shifth<0
    for i=0:-shifth-1
        x2a(1+i,1-min(0,shiftv):cv1-delta2-1) = x2(1,1+max(shiftv,0):cv1-delta2-1+shiftv);
        if shiftv<0
            for j=0:-shiftv-1
                x2a(1+i,1+j) = x2(1,1);
            end
        end
    end
end
if shiftv<0
    for i=0:-shiftv-1
        x2a(1-min(0,shifth):ch1-delta2-1,1+i) = x2(1+max(shifth,0):ch1-delta2-1+shifth,1);
    end
end


shifth = bestshifth(indh,1);
shiftv = bestshiftv(indh,1);
x2a(ch+delta2+1:h-max(0,shifth),1-min(0,shiftv):cv1-delta2-1) = x2(ch+shifth+delta2+1:h+min(shifth,0),1+max(shiftv,0):cv1-delta2-1+shiftv);
if shifth>0
    for i=0:shifth-1
        x2a(h-i,1-min(0,shiftv):cv1-delta2-1) = x2(h,1+max(shiftv,0):cv1-delta2-1+shiftv);
        if shiftv<0
            for j=0:-shiftv-1
                x2a(h-i,1+j) = x2(h,1);
            end
        end
    end
end
if shiftv<0
    for i=0:-shiftv-1
        x2a(ch+delta2+1:h-max(0,shifth),1+i) = x2(ch+shifth+delta2+1:h+min(shifth,0),1);
    end
end


shifth = bestshifth(1,indv);
shiftv = bestshiftv(1,indv);
x2a(1-min(0,shifth):ch1-delta2-1,cv+delta2+1:v-max(0,shiftv)) = x2(1+max(shifth,0):ch1-delta2-1+shifth,cv+shiftv+delta2+1:v+min(shiftv,0));
if shifth<0
    for i=0:-shifth-1
        x2a(1+i,cv+delta2+1:v-max(0,shiftv)) = x2(1,cv+shiftv+delta2+1:v+min(shiftv,0));
        if shiftv>0
            for j=0:-shiftv-1
                x2a(1+i,v-j) = x2(1,v);
            end
        end
    end
end
if shiftv>0
    for i=0:-shiftv-1
        x2a(1-min(0,shifth):ch1-delta2-1,v-i) = x2(1+max(shifth,0):ch1-delta2-1+shifth,v);
    end
end


shifth = bestshifth(indh,indv);
shiftv = bestshiftv(indh,indv);
x2a(ch+delta2+1:h-max(0,shifth),cv+delta2+1:v-max(0,shiftv)) = x2(ch+shifth+delta2+1:h+min(shifth,0),cv+shiftv+delta2+1:v+min(shiftv,0));
if shifth>0
    for i=0:shifth-1
        x2a(h-i,cv+delta2+1:v-max(0,shiftv)) = x2(h,cv+shiftv+delta2+1:v+min(shiftv,0));
        if shiftv>0
            for j=0:shiftv-1
                x2a(h-i,v-j) = x2(h,v);
            end
        end
    end
end
if shiftv>0
    for i=0:shiftv-1
        x2a(ch+delta2+1:h-max(0,shifth),v-i) = x2(ch+shifth+delta2+1:h+min(shifth,0),v);
    end
end
