%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 [P,enmin]=minSpanTreeTri_DP_SC_simi(X0,Y,lambda)
X=choose_refer_and_recirc(X0,'close');%%shift to set the reference point as the 1st point
Nx=size(X,1);Ny=size(Y,1);
%%%preserve mapping between X0 and X
dif=repmat(permute(X0,[1,3,2]),[1,Nx,1])-repmat(permute(X,[3,1,2]),[Nx,1,1]);
len=sqrt(sum(dif.^2,3));
orgcor=zeros(Nx);
for ii=1:Nx
    [tmp,ind]=min(len(ii,:));
    orgcor(ii,ind)=1;
end
%%%%%
scdist=SC_cost_fast_fun([(1:Nx)' ones(Nx,1)],X,Y,50,0);%%Ny x Ny x Nx, 50 angles for quantification
edges=minSpanTree(X(2:Nx,:));% Nx-1 x 2, edges of minimum spanning tree of X(2:end,:)
edges=edges+1;

dify=repmat(permute(Y,[1,3,2]),[1,Ny,1])-repmat(permute(Y,[3,1,2]),[Ny,1,1]);%%jj versus kk
leny=sqrt(sum(dify.^2,3));

Ypair=[repmat(permute(Y,[2,1]),[1,1,Ny]);repmat(permute(Y,[2,3,1]),[1,Ny,1])];%4 x Ny x Ny, point pair of Y
Ypaircon=reshape(Ypair,[4,Ny*Ny]);%4 x Ny*Ny, concatenation of Ypair

Jx=[permute(X,[2,3,1]),permute([-X(:,2) X(:,1)],[2,3,1]),repmat(eye(2),[1,1,Nx])];%[a -b;b a]*[xi;yi]+[d1;d2],2 x 4 x Nx similarity transform

%%nearest neighbour indices for each point in Y
Nn=15;%%number of nearest neighbours
neigh=nearest_neighbour(leny,Nn);%% Nn x Ny
% lenyneigh=zeros(Nn,Ny);%%distance for each point and its nearest neighbours
% for ii=1:Ny
%     lenyneigh(:,ii)=leny(neigh(:,ii),ii);
% end
%%%
distx=sqrt(sum((X(2:Nx,:)-ones(Nx-1,1)*X(1,:)).^2,2));
[tmp,ind]=min(distx);
root=ind+1;%root of the min spanning tree, range is [2, Nx]

pointer=zeros(Ny,Ny,Nx);%recorder 
[energy,pointer]=treeTri_DP_simi_faster2(root,edges,scdist,X,Jx,leny,Ypaircon,pointer,neigh,lambda);%,lenyneigh);
% [energy,pointer]=treeTri_DP_simi(root,edges,scdist,Jx,Ypaircon,pointer);
% lenxneig=sqrt(sumsqr(X(root,:)-X(1,:)));
% energy=energy+abs(leny-lenxneig);

%%%%%%%%%%recover correspondence
corr=zeros(Nx,1);
enmin=min(energy(:));
[indx,indy]=find(energy==enmin);
corr(1)=indy(1);%reference
corr(root)=indx(1);%root
pointer=permute(pointer(:,indy(1),:),[1,3,2]);%Ny x Nx
corr=tree_match_recover(root,edges,pointer,corr);
%%%%%%%%%%%%%%%%permutation matrix
P=zeros(Nx,Ny);
for ii=1:Nx
    P(ii,corr(ii))=1;
end
%%%%%%%%%%%%%%%%%%%%%%%%%
P=orgcor*P;
