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