function [Corr]=pointmatchfun(points,labels,maxpoints)
%maxpoints is the maximal number of points a single label can match to

% npoints=length(points);
% nlabels=length(labels);
% Corr = randperm(npoints);

npoints=length(points);
nlabels=length(labels);

ncor = npoints*nlabels;

A = calculateAffinityMatrix(points, labels);
x = ones(ncor,1);
x = x/sum(x);

% Old_Tar = x'*A*x;
Incre = 10;
% stop_th = 0.11;
stop_th = 0.000000001;
while Incre>stop_th
    Old_Tar = x'*A*x;
    x_new = x.*((A*x)/(x'*A*x));
    Incre = x_new'*A*x_new - Old_Tar;
    
    % decide early stop
%     is_stop = 0;
%     for i_p = 1:npoints
%         
%         x_p_start = (i_p-1) * nlabels + 1;
%         x_p_end   = x_p_start + nlabels - 1;
%         x_p = x(x_p_start:x_p_end);
%         if sum(x_p)==0
%             is_stop = 1;
%         end
%     end
%     if is_stop == 1
%         break;
%     end
    x = x_new;
end


Corr = zeros(1,npoints);
mark_label = zeros(1, nlabels);
% Corr(:) = nlabels + 1;
for i_p = 1:npoints
    [x_max x_idx] = max(x);
    
    if (x_max(1)>0)
        
        idx_p = floor((x_idx-1)/nlabels);
        idx_l = (x_idx-1)- idx_p*nlabels;
        
        Corr(idx_p+1) = idx_l+1;
        
        mark_label(idx_l+1) = mark_label(idx_l+1) + 1;
        
        if mark_label(idx_l+1) > maxpoints
            x(idx_l+1:nlabels:end) = 0;
        end
        x_p_start = idx_p * nlabels + 1;
        x_p_end   = x_p_start + nlabels - 1;
        x(x_p_start:x_p_end) = 0;
    end
end


%% internal functions
function M = calculateAffinityMatrix(P, Q)

num_point_p = size(P,1);
num_point_q = size(Q,1);

%sigmad =0.04;
sigmad = 0.5;
Scale  = 1;

num_vertex = num_point_p * num_point_q;

dist_point_p = zeros(num_point_p, num_point_p);

for i=1:num_point_p
    for j=1:num_point_p
        dist_point_p(i,j) =  dist2(P(i,:), P(j,:));
    end
end

dist_point_q = zeros(num_point_q, num_point_q);
for i=1:num_point_q
    for j=1:num_point_q
        dist_point_q(i,j) =  dist2(Q(i,:), Q(j,:));
    end
end

M = zeros(num_vertex, num_vertex);


for i_v = 0:num_vertex-2
    
    idx_p_i = floor(i_v / num_point_q);
    idx_q_i = i_v - idx_p_i * num_point_q;
    
    idx_p_i = idx_p_i + 1;
    idx_q_i = idx_q_i + 1;
    
    for j_v=i_v+1:1:num_vertex-1
        
        idx_p_j = floor(j_v / num_point_q);
        idx_q_j = j_v - idx_p_j * num_point_q;

        idx_p_j = idx_p_j + 1;
        idx_q_j = idx_q_j + 1;
        
        DistDiff = abs(dist_point_p(idx_p_i, idx_p_j) - Scale*dist_point_q(idx_q_i, idx_q_j));
        if(DistDiff>=3*sigmad)
             M(i_v+1, j_v+1) = 0;
        else
            M(i_v+1, j_v+1) = 4.5-0.5*DistDiff*DistDiff /sigmad/sigmad; 
            
        end
    end
end

for i_v = 1:num_vertex
    for j_v = 1:i_v-1
        M(i_v, j_v) = M(j_v, i_v);
    end
end


function dist_pi_pj = dist2(pi, pj)
    dist_pi_pj = sqrt(sum((pi-pj).^2));

