function [Corr]=pointmatchfun(points,labels, maxpoints)
    %maxpoints is the maximal number of points a single label can match to
    % jeff newell, jesse glass, jo catrambone
    %     
    npoints=length(points);
    nlabels=length(labels);
    n = npoints*nlabels;
    % build Affinity Matrix over point-label mapping combinations
    % the conversion from r (or symmetrically, column) is:
    %   labelidx = ceil(r/npoints), pointidx = add one modulo(r,numpoints)
    % or the other way: r = (labelidx - 1)*npoints + pointidx , a unique
    %       factorization.
    % and symmetrically for c 
    A = zeros(n);
    sigmaconst=1; 

    for r = 1:n
        for c = 1:n
            r_labelidx = ceil(r/npoints);
            r_pointidx = rem(r - 1,npoints) + 1;
            c_labelidx = ceil(c/npoints);
            c_pointidx = rem(c - 1,npoints) + 1;
            if r == c
                continue;
            else
                deltapntX = abs(points(r_pointidx,1) - points(c_pointidx,1));
                deltapntY = abs(points(r_pointidx,2) - points(c_pointidx,2));
                deltalblX = abs(labels(r_labelidx,1) - labels(c_labelidx,1));
                deltalblY = abs(labels(r_labelidx,2) - labels(c_labelidx,2));
                A(r,c) = exp( -abs( sqrt( (deltapntX)^2 + (deltapntY)^2) - sqrt((deltalblX)^2 + (deltalblY)^2)/sigmaconst));
            end
        end
    end

    %  next use A to find dense subgraph with the Replicator algorithm...
    %       we want a dense subgraph with npoints (which is >= nlabels)
    Delta = 0.00001; 
    X = 0.5*ones(n,1);
    for t = 0:10000 
        Xprev = X;
        v = (A*X)/(X'*A*X);
        X = bsxfun(@times,X,v); 
        if sum(abs(Xprev-X)) < Delta
            break;
        end
    end

    Corr = zeros(1, npoints);
    [~,sortIdx] = sort(X,'descend');
    sortIdx(npoints+1:end) = [];
    sortIdx = sort(sortIdx);
    for i=[1:npoints]
        r = sortIdx(i);
        pt = rem(r-1,npoints) + 1;
        lb = ceil(r/npoints);
        Corr(pt) = lb; 
    end
end

