OpenShot Library | libopenshot 0.2.7
Hungarian.cpp
Go to the documentation of this file.
1///////////////////////////////////////////////////////////////////////////////
2// Hungarian.cpp: Implementation file for Class HungarianAlgorithm.
3//
4// This is a C++ wrapper with slight modification of a hungarian algorithm implementation by Markus Buehren.
5// The original implementation is a few mex-functions for use in MATLAB, found here:
6// http://www.mathworks.com/matlabcentral/fileexchange/6543-functions-for-the-rectangular-assignment-problem
7//
8// Both this code and the orignal code are published under the BSD license.
9// by Cong Ma, 2016
10//
11
12#include "Hungarian.h"
13
14using namespace std;
15
18
19//********************************************************//
20// A single function wrapper for solving assignment problem.
21//********************************************************//
23 vector<vector<double>> &DistMatrix,
24 vector<int> &Assignment)
25{
26 unsigned int nRows = DistMatrix.size();
27 unsigned int nCols = DistMatrix[0].size();
28
29 double *distMatrixIn = new double[nRows * nCols];
30 int *assignment = new int[nRows];
31 double cost = 0.0;
32
33 // Fill in the distMatrixIn. Mind the index is "i + nRows * j".
34 // Here the cost matrix of size MxN is defined as a double precision array of N*M elements.
35 // In the solving functions matrices are seen to be saved MATLAB-internally in row-order.
36 // (i.e. the matrix [1 2; 3 4] will be stored as a vector [1 3 2 4], NOT [1 2 3 4]).
37 for (unsigned int i = 0; i < nRows; i++)
38 for (unsigned int j = 0; j < nCols; j++)
39 distMatrixIn[i + nRows * j] = DistMatrix[i][j];
40
41 // call solving function
42 assignmentoptimal(assignment, &cost, distMatrixIn, nRows, nCols);
43
44 Assignment.clear();
45 for (unsigned int r = 0; r < nRows; r++)
46 Assignment.push_back(assignment[r]);
47
48 delete[] distMatrixIn;
49 delete[] assignment;
50 return cost;
51}
52
53//********************************************************//
54// Solve optimal solution for assignment problem using Munkres algorithm, also known as Hungarian Algorithm.
55//********************************************************//
56void HungarianAlgorithm::assignmentoptimal(
57 int *assignment,
58 double *cost,
59 double *distMatrixIn,
60 int nOfRows,
61 int nOfColumns)
62{
63 double *distMatrix, *distMatrixTemp, *distMatrixEnd, *columnEnd, value, minValue;
64 bool *coveredColumns, *coveredRows, *starMatrix, *newStarMatrix, *primeMatrix;
65 int nOfElements, minDim, row, col;
66
67 /* initialization */
68 *cost = 0;
69 for (row = 0; row < nOfRows; row++)
70 assignment[row] = -1;
71
72 /* generate working copy of distance Matrix */
73 /* check if all matrix elements are positive */
74 nOfElements = nOfRows * nOfColumns;
75 distMatrix = (double *)malloc(nOfElements * sizeof(double));
76 distMatrixEnd = distMatrix + nOfElements;
77
78 for (row = 0; row < nOfElements; row++)
79 {
80 value = distMatrixIn[row];
81 if (value < 0)
82 cerr << "All matrix elements have to be non-negative." << endl;
83 distMatrix[row] = value;
84 }
85
86 /* memory allocation */
87 coveredColumns = (bool *)calloc(nOfColumns, sizeof(bool));
88 coveredRows = (bool *)calloc(nOfRows, sizeof(bool));
89 starMatrix = (bool *)calloc(nOfElements, sizeof(bool));
90 primeMatrix = (bool *)calloc(nOfElements, sizeof(bool));
91 newStarMatrix = (bool *)calloc(nOfElements, sizeof(bool)); /* used in step4 */
92
93 /* preliminary steps */
94 if (nOfRows <= nOfColumns)
95 {
96 minDim = nOfRows;
97
98 for (row = 0; row < nOfRows; row++)
99 {
100 /* find the smallest element in the row */
101 distMatrixTemp = distMatrix + row;
102 minValue = *distMatrixTemp;
103 distMatrixTemp += nOfRows;
104 while (distMatrixTemp < distMatrixEnd)
105 {
106 value = *distMatrixTemp;
107 if (value < minValue)
108 minValue = value;
109 distMatrixTemp += nOfRows;
110 }
111
112 /* subtract the smallest element from each element of the row */
113 distMatrixTemp = distMatrix + row;
114 while (distMatrixTemp < distMatrixEnd)
115 {
116 *distMatrixTemp -= minValue;
117 distMatrixTemp += nOfRows;
118 }
119 }
120
121 /* Steps 1 and 2a */
122 for (row = 0; row < nOfRows; row++)
123 for (col = 0; col < nOfColumns; col++)
124 if (fabs(distMatrix[row + nOfRows * col]) < DBL_EPSILON)
125 if (!coveredColumns[col])
126 {
127 starMatrix[row + nOfRows * col] = true;
128 coveredColumns[col] = true;
129 break;
130 }
131 }
132 else /* if(nOfRows > nOfColumns) */
133 {
134 minDim = nOfColumns;
135
136 for (col = 0; col < nOfColumns; col++)
137 {
138 /* find the smallest element in the column */
139 distMatrixTemp = distMatrix + nOfRows * col;
140 columnEnd = distMatrixTemp + nOfRows;
141
142 minValue = *distMatrixTemp++;
143 while (distMatrixTemp < columnEnd)
144 {
145 value = *distMatrixTemp++;
146 if (value < minValue)
147 minValue = value;
148 }
149
150 /* subtract the smallest element from each element of the column */
151 distMatrixTemp = distMatrix + nOfRows * col;
152 while (distMatrixTemp < columnEnd)
153 *distMatrixTemp++ -= minValue;
154 }
155
156 /* Steps 1 and 2a */
157 for (col = 0; col < nOfColumns; col++)
158 for (row = 0; row < nOfRows; row++)
159 if (fabs(distMatrix[row + nOfRows * col]) < DBL_EPSILON)
160 if (!coveredRows[row])
161 {
162 starMatrix[row + nOfRows * col] = true;
163 coveredColumns[col] = true;
164 coveredRows[row] = true;
165 break;
166 }
167 for (row = 0; row < nOfRows; row++)
168 coveredRows[row] = false;
169 }
170
171 /* move to step 2b */
172 step2b(assignment, distMatrix, starMatrix, newStarMatrix, primeMatrix, coveredColumns, coveredRows, nOfRows, nOfColumns, minDim);
173
174 /* compute cost and remove invalid assignments */
175 computeassignmentcost(assignment, cost, distMatrixIn, nOfRows);
176
177 /* free allocated memory */
178 free(distMatrix);
179 free(coveredColumns);
180 free(coveredRows);
181 free(starMatrix);
182 free(primeMatrix);
183 free(newStarMatrix);
184
185 return;
186}
187
188/********************************************************/
189void HungarianAlgorithm::buildassignmentvector(
190 int *assignment,
191 bool *starMatrix,
192 int nOfRows,
193 int nOfColumns)
194{
195 for (int row = 0; row < nOfRows; row++)
196 for (int col = 0; col < nOfColumns; col++)
197 if (starMatrix[row + nOfRows * col])
198 {
199#ifdef ONE_INDEXING
200 assignment[row] = col + 1; /* MATLAB-Indexing */
201#else
202 assignment[row] = col;
203#endif
204 break;
205 }
206}
207
208/********************************************************/
209void HungarianAlgorithm::computeassignmentcost(
210 int *assignment,
211 double *cost,
212 double *distMatrix,
213 int nOfRows)
214{
215 int row, col;
216
217 for (row = 0; row < nOfRows; row++)
218 {
219 col = assignment[row];
220 if (col >= 0)
221 *cost += distMatrix[row + nOfRows * col];
222 }
223}
224
225/********************************************************/
226void HungarianAlgorithm::step2a(
227 int *assignment,
228 double *distMatrix,
229 bool *starMatrix,
230 bool *newStarMatrix,
231 bool *primeMatrix,
232 bool *coveredColumns,
233 bool *coveredRows,
234 int nOfRows,
235 int nOfColumns,
236 int minDim)
237{
238 bool *starMatrixTemp, *columnEnd;
239 int col;
240
241 /* cover every column containing a starred zero */
242 for (col = 0; col < nOfColumns; col++)
243 {
244 starMatrixTemp = starMatrix + nOfRows * col;
245 columnEnd = starMatrixTemp + nOfRows;
246 while (starMatrixTemp < columnEnd)
247 {
248 if (*starMatrixTemp++)
249 {
250 coveredColumns[col] = true;
251 break;
252 }
253 }
254 }
255
256 /* move to step 3 */
257 step2b(assignment, distMatrix, starMatrix, newStarMatrix, primeMatrix, coveredColumns, coveredRows, nOfRows, nOfColumns, minDim);
258}
259
260/********************************************************/
261void HungarianAlgorithm::step2b(
262 int *assignment,
263 double *distMatrix,
264 bool *starMatrix,
265 bool *newStarMatrix,
266 bool *primeMatrix,
267 bool *coveredColumns,
268 bool *coveredRows,
269 int nOfRows,
270 int nOfColumns,
271 int minDim)
272{
273 int col, nOfCoveredColumns;
274
275 /* count covered columns */
276 nOfCoveredColumns = 0;
277 for (col = 0; col < nOfColumns; col++)
278 if (coveredColumns[col])
279 nOfCoveredColumns++;
280
281 if (nOfCoveredColumns == minDim)
282 {
283 /* algorithm finished */
284 buildassignmentvector(assignment, starMatrix, nOfRows, nOfColumns);
285 }
286 else
287 {
288 /* move to step 3 */
289 step3(assignment, distMatrix, starMatrix, newStarMatrix, primeMatrix, coveredColumns, coveredRows, nOfRows, nOfColumns, minDim);
290 }
291}
292
293/********************************************************/
294void HungarianAlgorithm::step3(
295 int *assignment,
296 double *distMatrix,
297 bool *starMatrix,
298 bool *newStarMatrix,
299 bool *primeMatrix,
300 bool *coveredColumns,
301 bool *coveredRows,
302 int nOfRows,
303 int nOfColumns,
304 int minDim)
305{
306 bool zerosFound;
307 int row, col, starCol;
308
309 zerosFound = true;
310 while (zerosFound)
311 {
312 zerosFound = false;
313 for (col = 0; col < nOfColumns; col++)
314 if (!coveredColumns[col])
315 for (row = 0; row < nOfRows; row++)
316 if ((!coveredRows[row]) && (fabs(distMatrix[row + nOfRows * col]) < DBL_EPSILON))
317 {
318 /* prime zero */
319 primeMatrix[row + nOfRows * col] = true;
320
321 /* find starred zero in current row */
322 for (starCol = 0; starCol < nOfColumns; starCol++)
323 if (starMatrix[row + nOfRows * starCol])
324 break;
325
326 if (starCol == nOfColumns) /* no starred zero found */
327 {
328 /* move to step 4 */
329 step4(assignment, distMatrix, starMatrix, newStarMatrix, primeMatrix, coveredColumns, coveredRows, nOfRows, nOfColumns, minDim, row, col);
330 return;
331 }
332 else
333 {
334 coveredRows[row] = true;
335 coveredColumns[starCol] = false;
336 zerosFound = true;
337 break;
338 }
339 }
340 }
341
342 /* move to step 5 */
343 step5(assignment, distMatrix, starMatrix, newStarMatrix, primeMatrix, coveredColumns, coveredRows, nOfRows, nOfColumns, minDim);
344}
345
346/********************************************************/
347void HungarianAlgorithm::step4(
348 int *assignment,
349 double *distMatrix,
350 bool *starMatrix,
351 bool *newStarMatrix,
352 bool *primeMatrix,
353 bool *coveredColumns,
354 bool *coveredRows,
355 int nOfRows,
356 int nOfColumns,
357 int minDim,
358 int row,
359 int col)
360{
361 int n, starRow, starCol, primeRow, primeCol;
362 int nOfElements = nOfRows * nOfColumns;
363
364 /* generate temporary copy of starMatrix */
365 for (n = 0; n < nOfElements; n++)
366 newStarMatrix[n] = starMatrix[n];
367
368 /* star current zero */
369 newStarMatrix[row + nOfRows * col] = true;
370
371 /* find starred zero in current column */
372 starCol = col;
373 for (starRow = 0; starRow < nOfRows; starRow++)
374 if (starMatrix[starRow + nOfRows * starCol])
375 break;
376
377 while (starRow < nOfRows)
378 {
379 /* unstar the starred zero */
380 newStarMatrix[starRow + nOfRows * starCol] = false;
381
382 /* find primed zero in current row */
383 primeRow = starRow;
384 for (primeCol = 0; primeCol < nOfColumns; primeCol++)
385 if (primeMatrix[primeRow + nOfRows * primeCol])
386 break;
387
388 /* star the primed zero */
389 newStarMatrix[primeRow + nOfRows * primeCol] = true;
390
391 /* find starred zero in current column */
392 starCol = primeCol;
393 for (starRow = 0; starRow < nOfRows; starRow++)
394 if (starMatrix[starRow + nOfRows * starCol])
395 break;
396 }
397
398 /* use temporary copy as new starMatrix */
399 /* delete all primes, uncover all rows */
400 for (n = 0; n < nOfElements; n++)
401 {
402 primeMatrix[n] = false;
403 starMatrix[n] = newStarMatrix[n];
404 }
405 for (n = 0; n < nOfRows; n++)
406 coveredRows[n] = false;
407
408 /* move to step 2a */
409 step2a(assignment, distMatrix, starMatrix, newStarMatrix, primeMatrix, coveredColumns, coveredRows, nOfRows, nOfColumns, minDim);
410}
411
412/********************************************************/
413void HungarianAlgorithm::step5(
414 int *assignment,
415 double *distMatrix,
416 bool *starMatrix,
417 bool *newStarMatrix,
418 bool *primeMatrix,
419 bool *coveredColumns,
420 bool *coveredRows,
421 int nOfRows,
422 int nOfColumns,
423 int minDim)
424{
425 double h, value;
426 int row, col;
427
428 /* find smallest uncovered element h */
429 h = DBL_MAX;
430 for (row = 0; row < nOfRows; row++)
431 if (!coveredRows[row])
432 for (col = 0; col < nOfColumns; col++)
433 if (!coveredColumns[col])
434 {
435 value = distMatrix[row + nOfRows * col];
436 if (value < h)
437 h = value;
438 }
439
440 /* add h to each covered row */
441 for (row = 0; row < nOfRows; row++)
442 if (coveredRows[row])
443 for (col = 0; col < nOfColumns; col++)
444 distMatrix[row + nOfRows * col] += h;
445
446 /* subtract h from each uncovered column */
447 for (col = 0; col < nOfColumns; col++)
448 if (!coveredColumns[col])
449 for (row = 0; row < nOfRows; row++)
450 distMatrix[row + nOfRows * col] -= h;
451
452 /* move to step 3 */
453 step3(assignment, distMatrix, starMatrix, newStarMatrix, primeMatrix, coveredColumns, coveredRows, nOfRows, nOfColumns, minDim);
454}
double Solve(std::vector< std::vector< double > > &DistMatrix, std::vector< int > &Assignment)
Definition: Hungarian.cpp:22