OpenShot Library | libopenshot 0.3.3
Loading...
Searching...
No Matches
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
18using 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//********************************************************//
60void 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/********************************************************/
193void 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/********************************************************/
213void 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/********************************************************/
230void 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/********************************************************/
265void 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/********************************************************/
298void 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/********************************************************/
351void 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/********************************************************/
417void 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}
double Solve(std::vector< std::vector< double > > &DistMatrix, std::vector< int > &Assignment)
Definition Hungarian.cpp:26