matthiasm@8: function [a,ass] = bipartiteMatchingIntProg(dst, nmatches)
matthiasm@8: % BIPARTITEMATCHINGINTPROG Use binary integer programming (linear objective) to solve for optimal linear assignment
matthiasm@8: % function a = bipartiteMatchingIntProg(dst)
matthiasm@8: % a(i) = best matching column for row i
matthiasm@8: % 
matthiasm@8: % This gives the same result as bipartiteMatchingHungarian.
matthiasm@8: %
matthiasm@8: % function a = bibpartiteMatchingIntProg(dst, nmatches)
matthiasm@8: % only matches the specified number (must be <= min(size(dst))).
matthiasm@8: % This can be used to allow outliers in both source and target.
matthiasm@8: %
matthiasm@8: % For details, see Marciel & Costeira, "A global solution to sparse correspondence
matthiasm@8: % problems", PAMI 25(2), 2003
matthiasm@8: 
matthiasm@8: if nargin < 2, nmatches = []; end
matthiasm@8: 
matthiasm@8: [p1 p2] = size(dst);
matthiasm@8: p1orig = p1; p2orig = p2;
matthiasm@8: dstorig = dst;
matthiasm@8: 
matthiasm@8: if isempty(nmatches) % no outliers allowed  (modulo size difference)
matthiasm@8:   % ensure matrix is square
matthiasm@8:   m = max(dst(:));
matthiasm@8:   if p1<p2
matthiasm@8:     dst = [dst; m*ones(p2-p1, p2)];
matthiasm@8:   elseif p1>p2
matthiasm@8:     dst = [dst  m*ones(p1, p1-p2)];
matthiasm@8:   end
matthiasm@8: end
matthiasm@8: [p1 p2] = size(dst);
matthiasm@8: 
matthiasm@8: 
matthiasm@8: c = dst(:); % vectorize cost matrix
matthiasm@8: 
matthiasm@8: % row-sum: ensure each column sums to 1
matthiasm@8: A2 = kron(eye(p2), ones(1,p1));
matthiasm@8: b2 = ones(p2,1);
matthiasm@8: 
matthiasm@8: % col-sum: ensure each row sums to 1
matthiasm@8: A3 = kron(ones(1,p2), eye(p1));
matthiasm@8: b3 = ones(p1,1);
matthiasm@8: 
matthiasm@8: if isempty(nmatches)
matthiasm@8:   % enforce doubly  stochastic
matthiasm@8:   A = [A2; A3];
matthiasm@8:   b = [b2; b3];
matthiasm@8:   Aineq = zeros(1, p1*p2);
matthiasm@8:   bineq = 0;
matthiasm@8: else
matthiasm@8:   nmatches = min([nmatches, p1, p2]);
matthiasm@8:   Aineq = [A2; A3];
matthiasm@8:   bineq = [b2; b3]; % row and col sums <= 1
matthiasm@8:   A = ones(1,p1*p2);
matthiasm@8:   b = nmatches; % total num matches = b (otherwise get degenerate soln)
matthiasm@8: end
matthiasm@8: 
matthiasm@8: 
matthiasm@8: ass = bintprog(c, Aineq, bineq, A, b);
matthiasm@8: ass = reshape(ass, p1, p2);
matthiasm@8: 
matthiasm@8: a = zeros(1, p1orig);
matthiasm@8: for i=1:p1orig
matthiasm@8:   ndx = find(ass(i,:)==1);
matthiasm@8:   if ~isempty(ndx) & (ndx <= p2orig)
matthiasm@8:     a(i) = ndx;
matthiasm@8:   end
matthiasm@8: end
matthiasm@8: 
matthiasm@8: