function compare( path1, path2 )
% compare( dir1, dir2 )
% dir1/2 : path of file for first/second frame
% 
% compares two frames.
% method : uses 'identify' to find unique tile neighborhoods
%   looks for identical neighborhoods between images
%   find transform between the two images
%   allows user to adjust matches between frames

% possible improvement: the way transform is handled. one would imagine
% that in the end, the last transform, all matches should be taken into
% account because all matches are correct. right now this isnt done,
% becuase it still considers each stage as equal to a first stage transform,
% where there could be many errors in the process. 


% use identify to find unique neighborhoods
figure(1)
title 'image1'
colormap gray
[com(1),s1] = identify( path1 );

figure(2)
title 'image2'
colormap gray
[com(2),s2] = identify( path2 );

if (s1 ~= s2)
    disp( 'large difference in number of tiles error, redoing')
    
    s = min(s1,s2);
    if s == s2
        figure(1)
        [com(1),s1] = identify(path1, s);
    else
        figure(2)
        [com(2),s2] = identify(path2, s);
    end
end

% load in necessary data
load( '-mat', path2, 'x','y', 'conn', 'N');
x2 = x;
y2 = y;
conn2 = conn;
N2 = N;
load( '-mat', path1, 'x','y', 'N', 'ymax', 'xmax', 'yave', 'conn');
image1 = zeros(ymax,xmax);

% u for unique
u = struct( 'ID', zeros(s1,0), 'n', [] , 'length', 0);
    
% find unique signatures from first image
for n = 1:com(1).length
    if com(1).count(n) == 1
        u.ID = [u.ID, com(1).ID(:,n)];
        u.n = [u.n, com(1).n(n)];
        u.length = u.length +1;
    end
end
    
% try to find a match for each unique tile signature to a tile in teh 2nd
% frame
match = zeros(N,1);
% match = [];
for n = 1:u.length
    % repeat current ID until is same size as common
    temp = repmat(u.ID(:,n), [1, com(2).length]);
    % compare current ID to each ID in common 2 list
    temp = xor( temp, com(2).ID );
    % sum up the results of each comparison
    % if sums to zero, perfect match
    temp = sum( temp );
    % gives column of perfect match if it exists, empty otherwise
    temp = find ( temp == 0 );
    
    if (numel(temp) == 1 & com(2).count(temp)==1)
        match(u.n(n)) = com(2).n(temp);
    end
end

% plot matches
figure(3)
graph( x,y,x2,y2 , match, N, image1, conn, conn2, N2 )

[match] = fix( x,y, x2,y2, match);

    
    
if length(find(match)) > 3
    transform = find_transform( match, x,y, x2,y2 , N);
else
    disp( ' too few matches' )
    return;
    % junk that probably doesnt work.
%     % compute vectors
%     vector = struct( 'x', [], 'y',[], 'mag' , [], 'dir', []);
%     for a = 1:l
%         % vector pointing from second frame to first frame item
%         vector(a).x = x(match(a,1)) - x2(match(a,2));
%         vector(a).y = y(match(a,1)) - y2(match(a,2));
%         vector(a).mag = (vector(a).x)^2 + (vector(a).y)^2;
%         vector(a).mag = sqrt(vector(a).mag);
%         % arc tan can only handle -90 to 90 degrees.
%         if vector(a).x > 0
%             % normal atan
%             vector(a).dir = atan( vector(a).y / vector(a).x );
%         else
%             % atan plus shift
%             vector(a).dir = atan( vector(a).y / vector(a).x ) + pi;
%         end
%     end
%     % average displacement vector
%     average = struct( 'x',sum([vector.x])/l , 'y',sum([vector.y])/l, ...
%         'mag',sum([vector.mag])/l , 'dir',sum([vector.dir])/l )
% 
%     % displacement vector (found by clustering magnitude and direction
%     % seperately)
%     displace = struct('x',[],'y',[],'mag',find_mag( vector ), 'dir', find_dir(vector));
%     displace.x = displace.mag * cos( displace.dir );
%     displace.y = displace.mag * sin( displace.dir );
%     displace 
% 
%     % kmeans clustering to find displacement
%     temp = [[vector.mag]' [vector.dir]'];
%     [temp,c] = kmeans(temp,2)

end

% try to transform image 2 to match figure 1
transform = reshape(transform, 2,3);
new = transform * [x2';y2';ones(1,numel(x2))];

% continue to correct errors and transform until user wants to stop
figure(5)
graph( x,y,new(1,:), new(2,:), match,N, image1, conn, conn2, N2 )
cont = input( 'fix again? [y/n]', 's');
while cont == 'y'
    % now allow for non-exact matches , by location / proximity
    N = numel(x);
    for a = 1:N
        n1 = findclosest( new(1,:), new(2,:), x(a),y(a));
        n2 = findclosest( x,y, new(1,n1), new(2,n1));
        if n2 == a 
            match(a) = n1;
        end
    end
    
    
    figure(4)
    graph( x,y,new(1,:), new(2,:), match,N, image1, conn, conn2, N2 )
    % correcting matches
    [match] = fix( x,y, new(1,:), new(2,:), match);
    
    transform = find_transform( match, x,y, new(1,:), new(2,:), N );
    % try to transform image 2 to match figure 1
    transform = reshape(transform, 2,3);
    new = transform * [new;ones(1,numel(new)/2)];
    
    
    figure(5)
    graph( x,y,new(1,:), new(2,:), match,N, image1, conn, conn2, N2 )
    
    cont = input( 'fix again? [y/n]', 's');
end

save( path1, 'match', '-append', '-mat');

clear
close all
%% ------------------------------------------------------------------------
function graph( x,y, x2,y2, m, N, im, conn, conn2, N2 )
% graphs connection graph of two frames on same image and the matches
% between them.


image(im)
colormap gray
hold on
% plot(x,y,'mx')
% plot(x2,y2,'c+')
hold off

if nargin > 7
    plotconn(x,y, conn, N, '.m-')
    plotconn(x2,y2, conn2, N2, '.c-')
end

hold on
clr = 'y-';
for a = 1:N
    if m(a) ~= 0
        plot( [x(a), x2(m(a))], [y(a), y2(m(a))], clr );
    end
end
hold off
%%-------------------------------------------------------------------------
function [match] = fix( x,y, xx2, yy2, match)
% allows user to modify matches


title 'LEFT click: create match, RIGHT click: remove match, CENTER click: exit/pause'

% correcting matches
disp( 'choose magenta then blue')
cont = 'y';
    while cont ~= 'q' 
        % get first point from user clicks
        [x1,y1,button] = ginput(1);
        
        if button == 2
            % pause or exit
            cont = input( '[p]ause [q]uit : ','s');
            if cont == 'p'
                pause;
            end
        else
            n1 = findclosest( x,y, floor(x1), floor(y1));
            hold on
            plot(x(n1), y(n1),'mo')

            % get second point from user clicks
            [x2,y2] = ginput(1);
            n2 = findclosest(xx2,yy2, floor(x2), floor(y2));
                plot( xx2(n2),yy2(n2),'co')
                
        % remove old matches
        if match(n1) ~= 0
            hold on 
            plot( [x(n1),xx2(match(n1))],[y(n1),yy2(match(n1))], 'k-');
            hold off
            match(n1) = 0;
        end
        rep = find(match == n2);
            for a = length(rep):-1:1
                r = rep(a);
                hold on
                plot( [x(r),xx2(match(r))],[y(r),yy2(match(r))], 'k-');
                hold off
                match(r) = 0;
            end
        end
                
        % add match
        if button == 1
            
            hold on
            plot( [x(n1), xx2(n2)], [y(n1), yy2(n2)], 'y-')
            hold off
            match(n1) = n2;
        end
  
    end
%%-------------------------------------------------------------------------
function transform = find_transform( m, x,y, x2,y2, l )
% transform = find_transform( vector, m, x,y, x2,y2 )
% m : stores coordinates of matches (previously 'match')

% threshold for accepting transforms
% should be low for small l and high for large l

T = [];
if l < 20
    % find all transforms
    % should be size [ l!/3!*(l-3)! , 6 ] order l^3
    for p1 = 1 : (l-2)
        if m(p1) ~= 0
        for p2 = (p1+1) : (l-1)
            if m(p2) ~= 0
            for p3 = (p2+1) : l
                if m(p3) ~= 0
                % points in first frame
                f1 = [ x(p1),     x(p2),     x(p3) ; ...
                       y(p1),     y(p2),     y(p3) ];
                % points in second frame
                f2 = [ x2(m(p1)), x2(m(p2)), x2(m(p3)) ; ...
                       y2(m(p1)), y2(m(p2)), y2(m(p3)) ; ...
                       1          , 1          , 1           ];
                % best fit transform matrix
                temp = f1 / f2 ;
                
                T = [ T ; temp(:)' ];
                end
            end 
            end
        end
        end
    end
else
    ll = length(find(m));
    if ll > 100
        lim = 16;
    elseif ll > 50
        lim = 8;
    else
        lim = 4;
    end
    
    option = find( m );
    oo = length(option);
    for b = 1 : oo
        p1 = option(b);
        if m(p1) ~= 0
            opt = option;
            opt(b) = [];
            % length of options
            o = oo-1;

            for a = 1:lim
                if o > 3 
                    if a > 1
                        r = randint(1,1,[1,o]);
                        p1 = opt( r );
                        opt(r) = [];
                        o = o - 1;
                    end

                    r = randint(1,1,[1,o]);
                    p2 = opt( r );
                    opt(r) = [];
                    o = o - 1;

                    r = randint(1,1,[1,o]);
                    p3 = opt( r );
                    opt(r) = [];
                    o = o - 1;       

                    % points in first frame
                    f1 = [ x(p1),  x(p2),  x(p3) ; ...
                           y(p1),  y(p2),  y(p3) ];
                    % points in second frame
                    f2 = [ x2(m(p1)), x2(m(p2)), x2(m(p3)) ; ...
                           y2(m(p1)), y2(m(p2)), y2(m(p3)) ; ...
                           1        , 1        , 1         ];
                    % best fit transform matrix
                    temp = f1 / f2 ;

                    T = [ T ; temp(:)' ];
                    end
            end
        end
    end
end
    

% ARBITRARY, the number 5 or 10 may work just as well.
% perhaps better in some cases (where most of the tiles should be taken
% into account, not just some)
group = ceil(sqrt(length(T)));
[ind, c] = kmeans ( T , group, 'replicates', 2);
% count the number of elements in each group
count = zeros(1,group);
for i = 1 : length(ind)
    count(ind(i)) = count(ind(i)) + 1;
end
[trash, i] = max (count);
% take cluster with most elements in it
transform = c(i,:);

%%-------------------------------------------------------------------------
function mag = find_mag( vector )
% mag = find_mag( vector )
% vector : array of vectors
% mag : magnitude found by closest cluster of magnitudes
% 
% try to find cluster of similar vectors --> this will be the desired
% magnitude (hopefully)

% sort magintude in ascending order
[temp,ind] = sort( [vector.mag], 'ascend');

% create matrix of the difference between adjacent magnitudes
temp2 = [temp,0] - [0,temp];
temp2(1) = [];
temp2(length(temp2)) =[];
len = length(temp2);

% find the section with the smallest differences bewteen mag,
sect = round(sqrt(len));
t = zeros(1,len-sect+1);
for ( i = 1:len)
    if i+sect - 1 <= len
        t(i) = sum(temp2(i:(i+sect-1)));
    end
end
[low,i] = min(t);

% if more than one section, look for larger section
while length(i)>1
    sect = ceil(sect*1.5);
    t = zeros(1,len);
    for ( i = 1:len)
        if i+sect - 1 <= len
            t(i) = sum(temp2(i:(i+sect-1)));
        end
    end
    [low,i] = min(t);
end    

here = i+floor(sect/2);
mag = vector(ind(here)).mag;
%%-------------------------------------------------------------------------
function dir = find_dir( vector )
% dir = find_dir( vector )
% vector : array of vectors
% dir : direction found by closest cluster of magnitudes (radians)
% 
% try to find cluster of similar vectors --> this will be the desired
% direction (hopefully)

% sort directions in ascending order
[temp,ind] = sort( [vector.dir], 'ascend');

% create matrix of the difference between adjacent directions
len = numel(temp);
temp2 = [temp,temp(1)] - [temp(len),temp];
temp2(len+1) =[];
% account for the fact that it loops around, and x = x+2pi
i = find( temp2 < 0 );
for h = 1: length(i)
    temp2(i(h)) = temp2(i(h)) + 2 * pi;
end

% find the section with the smallest differences bewteen dir.
sect = round(sqrt(len));
t = zeros(1,len);
for ( i = 1:len)
    if i+sect - 1 <= len
        t(i) = sum(temp2(i:(i+sect-1)));
    else
        t(i) = sum(temp2(i:len)) + sum(temp2(1:sect - (len-i+1)));
    end
end
[low,i] = min(t);

% if more than one section, look for larger section
while length(i)>1
    sect = ceil(sect*1.5);
    t = zeros(1,len);
    for ( i = 1:len)
        if i+sect - 1 <= len
            t(i) = sum(temp2(i:(i+sect-1)));
        else
            t(i) = sum(temp2(i:len)) + sum(temp2(1:sect - (len-i+1)));
        end
    end
    [low,i] = min(t);
end 

here = i+floor(sect/2);
dir = vector(ind(here)).dir;