%This software implements the point matching algorithm described in the paper:
%"Rotation Invariant Nonrigid Point Set Matching in Cluttered Scenes", IEEE transaction on Image Processing, 2012

%Author: Wei Lian (Changzhi University, Changzhi, 046011, Shanxi Province, P.R. China, Email: lianwei3 at gmail.com)

%If you use this software for research purposes, you should cite the aforementioned paper in any resulting publications.


function X=TSP_greedy(X)

% [X,Y]=generate_point_set('fishnoise',2); 


N=size(X,1);

dif=repmat(permute(X,[1,3,2]),[1,N,1])-repmat(permute(X,[3,1,2]),[N,1,1]);
dist=sqrt(sum(dif.^2,3));


for ii=1:N
    for jj=ii:N
        dist(ii,jj)=inf;
    end
end

edge=[];

while 1
    %%find the edge with mininal length
    min_dist=min(dist(:));
    [indx,indy]=find(dist==min_dist);
    indx=indx(1);
    indy=indy(1);
    %%set the edge to inf
    dist(indx,indy)=inf;
    %%check if the edge can be included to the TSP edge set
    xinedge=find(edge==indx);
    if size(xinedge)<2%%degree of vertex <2
        yinedge=find(edge==indy);
        if size(yinedge)<2
            if size(edge,1)==N-1%%this is the last edge
                edge=[edge;[indx,indy]];
                break;
            else
                %%check if there is a loop
                edge2=edge;
                indt=indx;
                while 1
                    [ind1,ind2]=find(edge2==indt);
                    if ind1
                        indt=edge2(ind1,3-ind2);%%get the index of the other vertex of the edge
                        edge2(ind1,:)=[inf,inf];
                    else
                        break;
                    end
                end
                if indt~=indy%%not a loop
                    edge=[edge;[indx, indy]];
                end
            end
        end
    end
end

%%convert edge expression to sequential list expression
Z=zeros(N,2);

ind1=1;
ind2=1;
ind=edge(1,1);
for ii=1:size(edge,1)     
    Z(ii,:)=X(ind,:);
    ind=edge(ind1,3-ind2);    
    edge(ind1,:)=[inf,inf];
    [ind1,ind2]=find(edge==ind);  
end


%%2-opt improvement
X=TSP_2opt_improvement(Z);%indispensable here
% X=Z;



