OpenShot Library | libopenshot  0.2.6
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 
14 using 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 //********************************************************//
56 void 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 /********************************************************/
189 void 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 /********************************************************/
209 void 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 /********************************************************/
226 void 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 /********************************************************/
261 void 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 /********************************************************/
294 void 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 /********************************************************/
347 void 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 /********************************************************/
413 void 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 }
STL namespace.
double Solve(std::vector< std::vector< double >> &DistMatrix, std::vector< int > &Assignment)
Definition: Hungarian.cpp:22