%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.


% addpath('../dynamicProg_shapeContext20091012');
% clear
% [X0,Y]=generate_point_set('fishnoise',2);%myplot(Y,'+');%return

function P=DPradiate_DPminSpanTreeSimi_orientedSC(X0,Y,alpha,lambda)
%This is a two stage algorithm.
%First, DPradiate is carried out to find the match of fan center.
%Then, assuming the match of fan center is known, DPboundary is
%carried out to find the matches of the rest points (i.e. boundary)
X=TSP_greedy(X0);%convert to closed contour
X=choose_refer_and_recirc(X,'close');%reorder to set the fan center 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));
cor=zeros(Nx);
for ii=1:Nx
    [tmp,ind]=min(len(ii,:));
    cor(ii,ind)=1;
end
%%%%%shape context distance
scdist=SC_cost_fast_fun([(1:Nx)' ones(Nx,1)],X,Y,50,0);%Ny x Ny x (Nx)
% alpha=.3;
%%difference vector for Y
dify=repmat(permute(Y,[1,3,2]),[1,Ny,1])-repmat(permute(Y,[3,1,2]),[Ny,1,1]);
leny=sqrt(sum(dify.^2,3));%Ny x Ny, length between every two points in Y
%%%%radiation shaped dynamic programming%%%stage one
% prev=zeros(Ny,Nx);
energy=zeros(1,Ny);
% alpha=1;
for ii=2:Nx
    lenx=sqrt(sumsqr(X(ii,:)-X(1,:)));%length of frame edge
    candidate=scdist(:,:,ii)+alpha*(leny-lenx).^2./(leny+lenx+eps);%abs(leny-lenx);
    candmin=min(candidate);
%     for jj=1:Ny
%         ind=find(candidate(:,jj)==candmin(jj));
%         prev(jj,ii)=ind(1);
%     end
    energy=energy+candmin;
end
%correspondence
corr=zeros(Nx,1);
% ind=find(energy==min(energy));
[tmp,ind]=min(energy);
corr(1)=ind;
% for ii=2:Nx
%     corr(ii)=prev(corr(1),ii);
% end
% %%%%%%%%%%%%%%%%permutation matrix
% P=zeros(Nx,Ny);
% for ii=1:Nx
%     P(ii,corr(ii))=1;
% end
% err=plot_intermediate(X,Y,P)
%%%%MST shaped dynamic programming%%%%stage two%%%%%
refer=corr(1);
edges=minSpanTree(X(2:Nx,:));% Nx-1 x 2, edges of minimum spanning tree of X(2:end,:)
edges=edges+1;
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
Ypair=[repmat(Y(refer,:)',[1,Ny]);Y'];%4 x Ny, point pair of Y, with Y(refer,:) known
pointer=zeros(Ny,Nx);%recorder 
root=2;
[energy,pointer]=treeTri_DP_simi_referenceKnown(root,edges,scdist,Jx,Ypair,pointer,refer,lambda);

%%%%%%%%%%recover correspondence
[tmp,indx]=min(energy);
corr(root)=indx;%root
% pointer=permute(pointer(:,refer,:),[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=cor*P;%%correct changing of correspondence caused by TPS process
% err=plot_intermediate(X0,Y,P)
