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