function [Corr]=pointmatchfun(points,labels, maxNum)

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

% Add labels when we have outliers 
ind = 0; 
max_value = max(max(max(points)), max(max(points))); 
if maxNum*size(labels,1)<size(points,1) 
    for i = 1:size(points,1)-maxNum*size(labels,1)
       labels = vertcat(labels, [1000*max_value 1000*max_value]); % One extra label for every outlier 
       ind = ind+1; 
    end    
end     

% We use Linear Programming as a Baseline predictor, the output is used to
% initialize the second method. 
% So first we will calculate the cost matrix 
% Cij is the "cost" of matching point i to job j 
C = [];
for i = 1:size(points,1)
Y = vertcat(points(i,:), labels);
D = squareform(pdist(Y));
D = D(1,2:end);
C = vertcat(C,D);  
end 
C = (reshape(C', [], 1))';
C =(C-min(min(C)))/(max(max(C))-min(min(C)));


u = [];
% Calculate the affinity matrix 
for i = 1:size(points,1)
    for j = 1:size(labels,1)
        u(end+1,:) = horzcat(points(i,:),  labels(j,:));
    end 
end 

v = u; 

for i = 1:size(u,1)
    for j = 1:size(v,1)
        if i==j 
            Affin(i,j) = 0;
        else  
            dp = sqrt(((u(i,1) - v(j,1))^2) + ((u(i,2) - v(j,2))^2)); 
            dl = sqrt(((u(i,3) - v(j,3))^2) + ((u(i,4) - v(j,4))^2)); 
            Affin(i,j) = exp(-(dp - dl)^2);
        end
    end 
end 

% Then you need to build the constraints:
A1 = zeros(size(points,1), size(points,1)*size(labels,1));
    for i = 1:size(points,1)
        A1(i, ((i-1)*size(labels,1)+1):(i)*size(labels,1) )  = 1;
    end 
b1 = ones(size(points,1), 1);

% 2. sumxij = 1 for all j : for all labels  
A2 = zeros(size(labels,1), size(points,1)*size(labels,1));
     for i = 1:size(labels,1) 
         A2(i, ((1:size(points,1))-1)*(size(labels,1))+i)  = 1;
     end 
b2 = ones(size(labels,1), 1)*maxNum;

% 3. Lower and Uper Bounds  
lb = zeros(size(labels,1)*size(points,1), 1);
ub = ones(size(labels,1)*size(points,1),1);

A = A2;
b = b2;
Aeq = A1;
Beq = b1;
 optimset('LargeScale','off','Simplex','on')
% Use simple linear programming to get the initial values 
[x0,fval0,exitflag0,output0]  = linprog(C',A,b,Aeq,Beq,lb,ub);
x0 = round(x0);
% exitflag0,output0, fval0

% This is the solution of the new method 
options = optimset('Algorithm','active-set', 'MaxIter',3000, 'PrecondBandWidth',inf);
[x,fval,exitflag,output]  = quadprog(-Affin,[],A,b,Aeq,Beq,lb,ub, x0, options); 
x = round(x);
% exitflag,output, fval
Corr = zeros(size(points,1), size(points,1)*size(labels,1));

% Prepare the output 
for i = 1:size(points,1)
Corr(i, ((i-1)*size(labels,1)+1):(i)*size(labels,1) )  = [1:size(labels,1)];
end 

Corr = (Corr*x)';
if ind>0
   for i = ind:1
    Corr(Corr(1,:) == size(points,1)-i+1) = 0; 
   end 
end 



