function [P,fup]=RPM_concave_point_matching(X,Y,err_dst,trans_type,alpha)
%X, Y: model and data point sets

%err_dst: average maximum allowed distance between a model point and its
%corresponding data point, which is later converted to tolerance error.

%trans_type: type of transformation, 'affine' for 2D affine or 'simi' for 2D similarity

%alpha: weight of regularization term. If affine transformation is
%employed, regularization is needed. If similarity transformation is
%employed, regularization is not needed, so alpha need not to be specified.


Nx=size(X,1);Ny=size(Y,1);

%contrary to the convention, P is converted to vector form by stacking its
%rows instead of its columns.

%%constraint for P
Arow=kron(eye(Nx),ones(1,Ny)); %% Nx x Nx*Ny, row equality constraints
Acol=kron(ones(1,Nx),eye(Ny)); %% Ny x Nx*Ny, column inequality constraints


switch trans_type
    case 'simi'
        %2D similarity transformation
        Jx=[X(:,1) -X(:,2) ones(Nx,1)*[1 0] X(:,2) X(:,1) ones(Nx,1)*[0 1]];%[x -y 1 0;y x 0 1], Jacobian
        Jx=reshape(Jx',[4,2*Nx])';%2*Nx x 4
        weight=(0)*diag([1 1 0 0]);%diagonal weight matrix
        trans0=[1 0 0 0]';%reference transformation

    case 'affine'
        %2D affine transformation
        Jx=[];
        for ii=1:size(X,1)
            Jx=[Jx;[kron(eye(size(X,2)),X(ii,:)) eye(size(X,2))]];
        end
%         alpha=1; %weight of the regularization term
        weight=alpha*diag([ones(1,size(X,2)^2) zeros(1,size(X,2))]);
        trans0=[reshape(eye(size(X,2)),[size(X,2)^2,1]); zeros(size(X,2),1)];
end

%convert RPM objective to quadratic form, here Asqr is the quadratic
%coefficient matrix, b is the linear coefficient vector
[Asqr,b]=RPM_objective_quadratic(Jx,Y,weight,trans0);

[Aq,Ar,ss]=qr(Asqr',0);
[U,V]=eigs(Ar*Ar',size(Jx,2));
U=Aq*U;
V=diag(V);


%branch and bound procedure
depth0=inf;%maximum allowed searching depth
err0=Nx*err_dst^2;%tolerance error
[P,fup]=normal_rectangular_branch_bound(Asqr,b,Acol,Arow,U,V,X,Y,depth0,err0,0);