%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 [energy,pointer]=treeTri_DP_simi_faster2(node,edges,scdist,X,Jx,leny,Ypaircon,pointer,neigh,lambda)         
%almost identical implementaion as treeTri_DP_simi_faster, the difference
%is at line 10, 31, 37 and 38
% lambda=1.;
Nn=size(neigh,1);
Ny=floor(sqrt(size(Ypaircon,2)));
lenxbase=sqrt(sumsqr(X(node,:)-X(1,:)));
energy=zeros(Ny,Ny);
%prediction step
invJnod=inv([Jx(:,:,node);Jx(:,:,1)]);%inv of Jacabian of point pair 1, node
transnod=reshape(invJnod*Ypaircon,[4,Ny,Ny]);%4 x Ny x Ny,similarity transformation between 1, ii and two points in Y
for ii=1:size(edges,1)
    if edges(ii,1)==node||edges(ii,2)==node
        if edges(ii,1)==node
            nex=edges(ii,2);
        else
            nex=edges(ii,1);
        end
        edges2=edges;
        edges2(ii,:)=[];%eliminate edge [node,nex]
        [leafen,pointer]=treeTri_DP_simi_faster2(nex,edges2,scdist,X,Jx,leny,Ypaircon,pointer,neigh,lambda);%leaf energy                                               
        %%%
        Nc=5;%%num of minimum points to the circle of some radius
        %%possible references for points in Y:
        %%points falling within the circular band with center some points in Y
        %%and radius equal to lenxbase and radial radius proportional to
        %%lenxbase
        bandrad=max(lenxbase/7,0.2);
        circband=nearest_neighbour_radius(abs(leny-lenxbase),bandrad,Nc);%%cell{Ny}(Nc), Nc is variable

        invJnex=inv([Jx(:,:,nex); Jx(:,:,1)]);%inv of Jacabian of point pair 1, nex
        transnex=reshape(invJnex*Ypaircon,[4,Ny,Ny]);
        
%         [leafen2, index]=DP_SC_fan_rotate_simi_speedup( leafen, transnod, transnex);
%         pointer(:,:,nex)=index;
        
        leafen2=inf*ones(Ny,Ny);
        for jj=1:Ny
            candidate=leafen(neigh(:,jj),circband{jj})+...
                lambda*squeeze(sum(abs(repmat(transnod(:,jj,circband{jj}),[1,Nn,1])-...
                transnex(:,neigh(:,jj),circband{jj})))); %cell{Ny}(Nn,Nc)                

            [sorted,index]=min(candidate);
            leafen2(jj,circband{jj})=sorted;
            for kk=1:size(circband{jj},1)
                pointer(jj,circband{jj}(kk),nex)=neigh(index(kk),jj);
            end            
        end
        energy=energy+leafen2;
    end
end
energy=energy+scdist(:,:,node);
