Code covered by the BSD License

# Functions for the rectangular assignment problem

### Markus Buehren (view profile)

14 Dec 2004 (Updated )

This package provides m- and mex-functions for solving the rectangular assignment problem.

assignmentoptimal(distMatrix)
```function [assignment, cost] = assignmentoptimal(distMatrix)
%ASSIGNMENTOPTIMAL    Compute optimal assignment by Munkres algorithm
%   ASSIGNMENTOPTIMAL(DISTMATRIX) computes the optimal assignment (minimum
%   overall costs) for the given rectangular distance or cost matrix, for
%   example the assignment of tracks (in rows) to observations (in
%   columns). The result is a column vector containing the assigned column
%   number in each row (or 0 if no assignment could be done).
%
%   [ASSIGNMENT, COST] = ASSIGNMENTOPTIMAL(DISTMATRIX) returns the
%   assignment vector and the overall cost.
%
%   The distance matrix may contain infinite values (forbidden
%   assignments). Internally, the infinite values are set to a very large
%   finite number, so that the Munkres algorithm itself works on
%   finite-number matrices. Before returning the assignment, all
%   assignments with infinite distance are deleted (i.e. set to zero).
%
%   A description of Munkres algorithm (also called Hungarian algorithm)
%   can easily be found on the web.
%
%   <a href="assignment.html">assignment.html</a>  <a href="http://www.mathworks.com/matlabcentral/fileexchange/6543">File Exchange</a>  <a href="https://www.paypal.com/cgi-bin/webscr?cmd=_s-xclick&hosted_button_id=EVW2A4G2HBVAU">Donate via PayPal</a>
%
%   Markus Buehren

% save original distMatrix for cost computation
originalDistMatrix = distMatrix;

% check for negative elements
if any(distMatrix(:) < 0)
error('All matrix elements have to be non-negative.');
end

% get matrix dimensions
[nOfRows, nOfColumns] = size(distMatrix);

% check for infinite values
finiteIndex   = isfinite(distMatrix);
infiniteIndex = find(~finiteIndex);
if ~isempty(infiniteIndex)
% set infinite values to large finite value
maxFiniteValue = max(max(distMatrix(finiteIndex)));
if maxFiniteValue > 0
infValue = abs(10 * maxFiniteValue * nOfRows * nOfColumns);
else
infValue = 10;
end
if isempty(infValue)
% all elements are infinite
assignment = zeros(nOfRows, 1);
cost       = 0;
return
end
distMatrix(infiniteIndex) = infValue;
end

% memory allocation
coveredColumns = zeros(1,       nOfColumns);
coveredRows    = zeros(nOfRows, 1);
starMatrix     = zeros(nOfRows, nOfColumns);
primeMatrix    = zeros(nOfRows, nOfColumns);

% preliminary steps
if nOfRows <= nOfColumns
minDim = nOfRows;

% find the smallest element of each row
minVector = min(distMatrix, [], 2);

% subtract the smallest element of each row from the row
distMatrix = distMatrix - repmat(minVector, 1, nOfColumns);

% Steps 1 and 2
for row = 1:nOfRows
for col = find(distMatrix(row,:)==0)
if ~coveredColumns(col)%~any(starMatrix(:,col))
starMatrix(row, col) = 1;
coveredColumns(col)  = 1;
break
end
end
end

else % nOfRows > nOfColumns
minDim = nOfColumns;

% find the smallest element of each column
minVector = min(distMatrix);

% subtract the smallest element of each column from the column
distMatrix = distMatrix - repmat(minVector, nOfRows, 1);

% Steps 1 and 2
for col = 1:nOfColumns
for row = find(distMatrix(:,col)==0)'
if ~coveredRows(row)
starMatrix(row, col) = 1;
coveredColumns(col)  = 1;
coveredRows(row)     = 1;
break
end
end
end
coveredRows(:) = 0; % was used auxiliary above
end

if sum(coveredColumns) == minDim
% algorithm finished
assignment = buildassignmentvector(starMatrix);
else
% move to step 3
[assignment, distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows] = ...
step3(distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows, minDim); %#ok
end

% compute cost and remove invalid assignments
[assignment, cost] = computeassignmentcost(assignment, originalDistMatrix, nOfRows);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function assignment = buildassignmentvector(starMatrix)

[maxValue, assignment] = max(starMatrix, [], 2);
assignment(maxValue == 0) = 0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [assignment, cost] = computeassignmentcost(assignment, distMatrix, nOfRows)

rowIndex   = find(assignment);
costVector = distMatrix(rowIndex + nOfRows * (assignment(rowIndex)-1));
finiteIndex = isfinite(costVector);
cost = sum(costVector(finiteIndex));
assignment(rowIndex(~finiteIndex)) = 0;

% Step 2: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [assignment, distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows] = ...
step2(distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows, minDim)

% cover every column containing a starred zero
maxValue = max(starMatrix);
coveredColumns(maxValue == 1) = 1;

if sum(coveredColumns) == minDim
% algorithm finished
assignment = buildassignmentvector(starMatrix);
else
% move to step 3
[assignment, distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows] = ...
step3(distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows, minDim);
end

% Step 3: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [assignment, distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows] = ...
step3(distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows, minDim)

zerosFound = 1;
while zerosFound

zerosFound = 0;
for col = find(~coveredColumns)
for row = find(~coveredRows')
if distMatrix(row,col) == 0

primeMatrix(row, col) = 1;
starCol = find(starMatrix(row,:));
if isempty(starCol)
% move to step 4
[assignment, distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows] = ...
step4(distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows, row, col, minDim);
return
else
coveredRows(row)        = 1;
coveredColumns(starCol) = 0;
zerosFound              = 1;
break % go on in next column
end
end
end
end
end

% move to step 5
[assignment, distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows] = ...
step5(distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows, minDim);

% Step 4: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [assignment, distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows] = ...
step4(distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows, row, col, minDim)

newStarMatrix          = starMatrix;
newStarMatrix(row,col) = 1;

starCol = col;
starRow = find(starMatrix(:, starCol));

while ~isempty(starRow)

% unstar the starred zero
newStarMatrix(starRow, starCol) = 0;

% find primed zero in row
primeRow = starRow;
primeCol = find(primeMatrix(primeRow, :));

% star the primed zero
newStarMatrix(primeRow, primeCol) = 1;

% find starred zero in column
starCol = primeCol;
starRow = find(starMatrix(:, starCol));

end
starMatrix = newStarMatrix;

primeMatrix(:) = 0;
coveredRows(:) = 0;

% move to step 2
[assignment, distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows] = ...
step2(distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows, minDim);

% Step 5: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [assignment, distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows] = ...
step5(distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows, minDim)

% find smallest uncovered element
uncoveredRowsIndex    = find(~coveredRows');
uncoveredColumnsIndex = find(~coveredColumns);
[s, index1] = min(distMatrix(uncoveredRowsIndex,uncoveredColumnsIndex));
[s, index2] = min(s); %#ok
h = distMatrix(uncoveredRowsIndex(index1(index2)), uncoveredColumnsIndex(index2));

% add h to each covered row
index = find(coveredRows);
distMatrix(index, :) = distMatrix(index, :) + h;

% subtract h from each uncovered column
distMatrix(:, uncoveredColumnsIndex) = distMatrix(:, uncoveredColumnsIndex) - h;

% move to step 3
[assignment, distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows] = ...
step3(distMatrix, starMatrix, primeMatrix, coveredColumns, coveredRows, minDim);

```