···11+///////////////////////////////////////////////////////////////////////////////
22+// Hungarian.h: Header file for Class HungarianAlgorithm.
33+//
44+// This is a C++ wrapper with slight modification of a hungarian algorithm implementation by Markus Buehren.
55+// The original implementation is a few mex-functions for use in MATLAB, found here:
66+// http://www.mathworks.com/matlabcentral/fileexchange/6543-functions-for-the-rectangular-assignment-problem
77+//
88+// Both this code and the orignal code are published under the BSD license.
99+// by Cong Ma, 2016
1010+//
1111+1212+#ifndef HUNGARIAN_H
1313+#define HUNGARIAN_H
1414+1515+#include <iostream>
1616+#include <vector>
1717+1818+using namespace std;
1919+2020+2121+class HungarianAlgorithm
2222+{
2323+public:
2424+ HungarianAlgorithm();
2525+ ~HungarianAlgorithm();
2626+ double Solve(vector <vector<double> >& DistMatrix, vector<int>& Assignment);
2727+2828+private:
2929+ void assignmentoptimal(int *assignment, double *cost, double *distMatrix, int nOfRows, int nOfColumns);
3030+ void buildassignmentvector(int *assignment, bool *starMatrix, int nOfRows, int nOfColumns);
3131+ void computeassignmentcost(int *assignment, double *cost, double *distMatrix, int nOfRows);
3232+ void step2a(int *assignment, double *distMatrix, bool *starMatrix, bool *newStarMatrix, bool *primeMatrix, bool *coveredColumns, bool *coveredRows, int nOfRows, int nOfColumns, int minDim);
3333+ void step2b(int *assignment, double *distMatrix, bool *starMatrix, bool *newStarMatrix, bool *primeMatrix, bool *coveredColumns, bool *coveredRows, int nOfRows, int nOfColumns, int minDim);
3434+ void step3(int *assignment, double *distMatrix, bool *starMatrix, bool *newStarMatrix, bool *primeMatrix, bool *coveredColumns, bool *coveredRows, int nOfRows, int nOfColumns, int minDim);
3535+ void step4(int *assignment, double *distMatrix, bool *starMatrix, bool *newStarMatrix, bool *primeMatrix, bool *coveredColumns, bool *coveredRows, int nOfRows, int nOfColumns, int minDim, int row, int col);
3636+ void step5(int *assignment, double *distMatrix, bool *starMatrix, bool *newStarMatrix, bool *primeMatrix, bool *coveredColumns, bool *coveredRows, int nOfRows, int nOfColumns, int minDim);
3737+};
3838+3939+4040+#endif///////////////////////////////////////////////////////////////////////////////
4141+// Hungarian.cpp: Implementation file for Class HungarianAlgorithm.
4242+//
4343+// This is a C++ wrapper with slight modification of a hungarian algorithm implementation by Markus Buehren.
4444+// The original implementation is a few mex-functions for use in MATLAB, found here:
4545+// http://www.mathworks.com/matlabcentral/fileexchange/6543-functions-for-the-rectangular-assignment-problem
4646+//
4747+// Both this code and the orignal code are published under the BSD license.
4848+// by Cong Ma, 2016
4949+//
5050+5151+#include <stdlib.h>
5252+#include <cfloat> // for DBL_MAX
5353+#include <cmath> // for fabs()
5454+//#include "Hungarian.hpp"
5555+5656+5757+HungarianAlgorithm::HungarianAlgorithm(){}
5858+HungarianAlgorithm::~HungarianAlgorithm(){}
5959+6060+6161+//********************************************************//
6262+// A single function wrapper for solving assignment problem.
6363+//********************************************************//
6464+double HungarianAlgorithm::Solve(vector <vector<double> >& DistMatrix, vector<int>& Assignment)
6565+{
6666+ unsigned int nRows = DistMatrix.size();
6767+ unsigned int nCols = DistMatrix[0].size();
6868+6969+ double *distMatrixIn = new double[nRows * nCols];
7070+ int *assignment = new int[nRows];
7171+ double cost = 0.0;
7272+7373+ // Fill in the distMatrixIn. Mind the index is "i + nRows * j".
7474+ // Here the cost matrix of size MxN is defined as a double precision array of N*M elements.
7575+ // In the solving functions matrices are seen to be saved MATLAB-internally in row-order.
7676+ // (i.e. the matrix [1 2; 3 4] will be stored as a vector [1 3 2 4], NOT [1 2 3 4]).
7777+ for (unsigned int i = 0; i < nRows; i++)
7878+ for (unsigned int j = 0; j < nCols; j++)
7979+ distMatrixIn[i + nRows * j] = DistMatrix[i][j];
8080+8181+ // call solving function
8282+ assignmentoptimal(assignment, &cost, distMatrixIn, nRows, nCols);
8383+8484+ Assignment.clear();
8585+ for (unsigned int r = 0; r < nRows; r++)
8686+ Assignment.push_back(assignment[r]);
8787+8888+ delete[] distMatrixIn;
8989+ delete[] assignment;
9090+ return cost;
9191+}
9292+9393+9494+//********************************************************//
9595+// Solve optimal solution for assignment problem using Munkres algorithm, also known as Hungarian Algorithm.
9696+//********************************************************//
9797+void HungarianAlgorithm::assignmentoptimal(int *assignment, double *cost, double *distMatrixIn, int nOfRows, int nOfColumns)
9898+{
9999+ double *distMatrix, *distMatrixTemp, *distMatrixEnd, *columnEnd, value, minValue;
100100+ bool *coveredColumns, *coveredRows, *starMatrix, *newStarMatrix, *primeMatrix;
101101+ int nOfElements, minDim, row, col;
102102+103103+ /* initialization */
104104+ *cost = 0;
105105+ for (row = 0; row<nOfRows; row++)
106106+ assignment[row] = -1;
107107+108108+ /* generate working copy of distance Matrix */
109109+ /* check if all matrix elements are positive */
110110+ nOfElements = nOfRows * nOfColumns;
111111+ distMatrix = (double *)malloc(nOfElements * sizeof(double));
112112+ distMatrixEnd = distMatrix + nOfElements;
113113+114114+ for (row = 0; row<nOfElements; row++)
115115+ {
116116+ value = distMatrixIn[row];
117117+ if (value < 0)
118118+ cerr << "All matrix elements have to be non-negative." << endl;
119119+ distMatrix[row] = value;
120120+ }
121121+122122+123123+ /* memory allocation */
124124+ coveredColumns = (bool *)calloc(nOfColumns, sizeof(bool));
125125+ coveredRows = (bool *)calloc(nOfRows, sizeof(bool));
126126+ starMatrix = (bool *)calloc(nOfElements, sizeof(bool));
127127+ primeMatrix = (bool *)calloc(nOfElements, sizeof(bool));
128128+ newStarMatrix = (bool *)calloc(nOfElements, sizeof(bool)); /* used in step4 */
129129+130130+ /* preliminary steps */
131131+ if (nOfRows <= nOfColumns)
132132+ {
133133+ minDim = nOfRows;
134134+135135+ for (row = 0; row<nOfRows; row++)
136136+ {
137137+ /* find the smallest element in the row */
138138+ distMatrixTemp = distMatrix + row;
139139+ minValue = *distMatrixTemp;
140140+ distMatrixTemp += nOfRows;
141141+ while (distMatrixTemp < distMatrixEnd)
142142+ {
143143+ value = *distMatrixTemp;
144144+ if (value < minValue)
145145+ minValue = value;
146146+ distMatrixTemp += nOfRows;
147147+ }
148148+149149+ /* subtract the smallest element from each element of the row */
150150+ distMatrixTemp = distMatrix + row;
151151+ while (distMatrixTemp < distMatrixEnd)
152152+ {
153153+ *distMatrixTemp -= minValue;
154154+ distMatrixTemp += nOfRows;
155155+ }
156156+ }
157157+158158+ /* Steps 1 and 2a */
159159+ for (row = 0; row<nOfRows; row++)
160160+ for (col = 0; col<nOfColumns; col++)
161161+ if (fabs(distMatrix[row + nOfRows*col]) < DBL_EPSILON)
162162+ if (!coveredColumns[col])
163163+ {
164164+ starMatrix[row + nOfRows*col] = true;
165165+ coveredColumns[col] = true;
166166+ break;
167167+ }
168168+ }
169169+ else /* if(nOfRows > nOfColumns) */
170170+ {
171171+ minDim = nOfColumns;
172172+173173+ for (col = 0; col<nOfColumns; col++)
174174+ {
175175+ /* find the smallest element in the column */
176176+ distMatrixTemp = distMatrix + nOfRows*col;
177177+ columnEnd = distMatrixTemp + nOfRows;
178178+179179+ minValue = *distMatrixTemp++;
180180+ while (distMatrixTemp < columnEnd)
181181+ {
182182+ value = *distMatrixTemp++;
183183+ if (value < minValue)
184184+ minValue = value;
185185+ }
186186+187187+ /* subtract the smallest element from each element of the column */
188188+ distMatrixTemp = distMatrix + nOfRows*col;
189189+ while (distMatrixTemp < columnEnd)
190190+ *distMatrixTemp++ -= minValue;
191191+ }
192192+193193+ /* Steps 1 and 2a */
194194+ for (col = 0; col<nOfColumns; col++)
195195+ for (row = 0; row<nOfRows; row++)
196196+ if (fabs(distMatrix[row + nOfRows*col]) < DBL_EPSILON)
197197+ if (!coveredRows[row])
198198+ {
199199+ starMatrix[row + nOfRows*col] = true;
200200+ coveredColumns[col] = true;
201201+ coveredRows[row] = true;
202202+ break;
203203+ }
204204+ for (row = 0; row<nOfRows; row++)
205205+ coveredRows[row] = false;
206206+207207+ }
208208+209209+ /* move to step 2b */
210210+ step2b(assignment, distMatrix, starMatrix, newStarMatrix, primeMatrix, coveredColumns, coveredRows, nOfRows, nOfColumns, minDim);
211211+212212+ /* compute cost and remove invalid assignments */
213213+ computeassignmentcost(assignment, cost, distMatrixIn, nOfRows);
214214+215215+ /* free allocated memory */
216216+ free(distMatrix);
217217+ free(coveredColumns);
218218+ free(coveredRows);
219219+ free(starMatrix);
220220+ free(primeMatrix);
221221+ free(newStarMatrix);
222222+223223+ return;
224224+}
225225+226226+/********************************************************/
227227+void HungarianAlgorithm::buildassignmentvector(int *assignment, bool *starMatrix, int nOfRows, int nOfColumns)
228228+{
229229+ int row, col;
230230+231231+ for (row = 0; row<nOfRows; row++)
232232+ for (col = 0; col<nOfColumns; col++)
233233+ if (starMatrix[row + nOfRows*col])
234234+ {
235235+#ifdef ONE_INDEXING
236236+ assignment[row] = col + 1; /* MATLAB-Indexing */
237237+#else
238238+ assignment[row] = col;
239239+#endif
240240+ break;
241241+ }
242242+}
243243+244244+/********************************************************/
245245+void HungarianAlgorithm::computeassignmentcost(int *assignment, double *cost, double *distMatrix, int nOfRows)
246246+{
247247+ int row, col;
248248+249249+ for (row = 0; row<nOfRows; row++)
250250+ {
251251+ col = assignment[row];
252252+ if (col >= 0)
253253+ *cost += distMatrix[row + nOfRows*col];
254254+ }
255255+}
256256+257257+/********************************************************/
258258+void HungarianAlgorithm::step2a(int *assignment, double *distMatrix, bool *starMatrix, bool *newStarMatrix, bool *primeMatrix, bool *coveredColumns, bool *coveredRows, int nOfRows, int nOfColumns, int minDim)
259259+{
260260+ bool *starMatrixTemp, *columnEnd;
261261+ int col;
262262+263263+ /* cover every column containing a starred zero */
264264+ for (col = 0; col<nOfColumns; col++)
265265+ {
266266+ starMatrixTemp = starMatrix + nOfRows*col;
267267+ columnEnd = starMatrixTemp + nOfRows;
268268+ while (starMatrixTemp < columnEnd){
269269+ if (*starMatrixTemp++)
270270+ {
271271+ coveredColumns[col] = true;
272272+ break;
273273+ }
274274+ }
275275+ }
276276+277277+ /* move to step 3 */
278278+ step2b(assignment, distMatrix, starMatrix, newStarMatrix, primeMatrix, coveredColumns, coveredRows, nOfRows, nOfColumns, minDim);
279279+}
280280+281281+/********************************************************/
282282+void HungarianAlgorithm::step2b(int *assignment, double *distMatrix, bool *starMatrix, bool *newStarMatrix, bool *primeMatrix, bool *coveredColumns, bool *coveredRows, int nOfRows, int nOfColumns, int minDim)
283283+{
284284+ int col, nOfCoveredColumns;
285285+286286+ /* count covered columns */
287287+ nOfCoveredColumns = 0;
288288+ for (col = 0; col<nOfColumns; col++)
289289+ if (coveredColumns[col])
290290+ nOfCoveredColumns++;
291291+292292+ if (nOfCoveredColumns == minDim)
293293+ {
294294+ /* algorithm finished */
295295+ buildassignmentvector(assignment, starMatrix, nOfRows, nOfColumns);
296296+ }
297297+ else
298298+ {
299299+ /* move to step 3 */
300300+ step3(assignment, distMatrix, starMatrix, newStarMatrix, primeMatrix, coveredColumns, coveredRows, nOfRows, nOfColumns, minDim);
301301+ }
302302+303303+}
304304+305305+/********************************************************/
306306+void HungarianAlgorithm::step3(int *assignment, double *distMatrix, bool *starMatrix, bool *newStarMatrix, bool *primeMatrix, bool *coveredColumns, bool *coveredRows, int nOfRows, int nOfColumns, int minDim)
307307+{
308308+ bool zerosFound;
309309+ int row, col, starCol;
310310+311311+ zerosFound = true;
312312+ while (zerosFound)
313313+ {
314314+ zerosFound = false;
315315+ for (col = 0; col<nOfColumns; col++)
316316+ if (!coveredColumns[col])
317317+ for (row = 0; row<nOfRows; row++)
318318+ if ((!coveredRows[row]) && (fabs(distMatrix[row + nOfRows*col]) < DBL_EPSILON))
319319+ {
320320+ /* prime zero */
321321+ primeMatrix[row + nOfRows*col] = true;
322322+323323+ /* find starred zero in current row */
324324+ for (starCol = 0; starCol<nOfColumns; starCol++)
325325+ if (starMatrix[row + nOfRows*starCol])
326326+ break;
327327+328328+ if (starCol == nOfColumns) /* no starred zero found */
329329+ {
330330+ /* move to step 4 */
331331+ step4(assignment, distMatrix, starMatrix, newStarMatrix, primeMatrix, coveredColumns, coveredRows, nOfRows, nOfColumns, minDim, row, col);
332332+ return;
333333+ }
334334+ else
335335+ {
336336+ coveredRows[row] = true;
337337+ coveredColumns[starCol] = false;
338338+ zerosFound = true;
339339+ break;
340340+ }
341341+ }
342342+ }
343343+344344+ /* move to step 5 */
345345+ step5(assignment, distMatrix, starMatrix, newStarMatrix, primeMatrix, coveredColumns, coveredRows, nOfRows, nOfColumns, minDim);
346346+}
347347+348348+/********************************************************/
349349+void HungarianAlgorithm::step4(int *assignment, double *distMatrix, bool *starMatrix, bool *newStarMatrix, bool *primeMatrix, bool *coveredColumns, bool *coveredRows, int nOfRows, int nOfColumns, int minDim, int row, int col)
350350+{
351351+ int n, starRow, starCol, primeRow, primeCol;
352352+ int nOfElements = nOfRows*nOfColumns;
353353+354354+ /* generate temporary copy of starMatrix */
355355+ for (n = 0; n<nOfElements; n++)
356356+ newStarMatrix[n] = starMatrix[n];
357357+358358+ /* star current zero */
359359+ newStarMatrix[row + nOfRows*col] = true;
360360+361361+ /* find starred zero in current column */
362362+ starCol = col;
363363+ for (starRow = 0; starRow<nOfRows; starRow++)
364364+ if (starMatrix[starRow + nOfRows*starCol])
365365+ break;
366366+367367+ while (starRow<nOfRows)
368368+ {
369369+ /* unstar the starred zero */
370370+ newStarMatrix[starRow + nOfRows*starCol] = false;
371371+372372+ /* find primed zero in current row */
373373+ primeRow = starRow;
374374+ for (primeCol = 0; primeCol<nOfColumns; primeCol++)
375375+ if (primeMatrix[primeRow + nOfRows*primeCol])
376376+ break;
377377+378378+ /* star the primed zero */
379379+ newStarMatrix[primeRow + nOfRows*primeCol] = true;
380380+381381+ /* find starred zero in current column */
382382+ starCol = primeCol;
383383+ for (starRow = 0; starRow<nOfRows; starRow++)
384384+ if (starMatrix[starRow + nOfRows*starCol])
385385+ break;
386386+ }
387387+388388+ /* use temporary copy as new starMatrix */
389389+ /* delete all primes, uncover all rows */
390390+ for (n = 0; n<nOfElements; n++)
391391+ {
392392+ primeMatrix[n] = false;
393393+ starMatrix[n] = newStarMatrix[n];
394394+ }
395395+ for (n = 0; n<nOfRows; n++)
396396+ coveredRows[n] = false;
397397+398398+ /* move to step 2a */
399399+ step2a(assignment, distMatrix, starMatrix, newStarMatrix, primeMatrix, coveredColumns, coveredRows, nOfRows, nOfColumns, minDim);
400400+}
401401+402402+/********************************************************/
403403+void HungarianAlgorithm::step5(int *assignment, double *distMatrix, bool *starMatrix, bool *newStarMatrix, bool *primeMatrix, bool *coveredColumns, bool *coveredRows, int nOfRows, int nOfColumns, int minDim)
404404+{
405405+ double h, value;
406406+ int row, col;
407407+408408+ /* find smallest uncovered element h */
409409+ h = DBL_MAX;
410410+ for (row = 0; row<nOfRows; row++)
411411+ if (!coveredRows[row])
412412+ for (col = 0; col<nOfColumns; col++)
413413+ if (!coveredColumns[col])
414414+ {
415415+ value = distMatrix[row + nOfRows*col];
416416+ if (value < h)
417417+ h = value;
418418+ }
419419+420420+ /* add h to each covered row */
421421+ for (row = 0; row<nOfRows; row++)
422422+ if (coveredRows[row])
423423+ for (col = 0; col<nOfColumns; col++)
424424+ distMatrix[row + nOfRows*col] += h;
425425+426426+ /* subtract h from each uncovered column */
427427+ for (col = 0; col<nOfColumns; col++)
428428+ if (!coveredColumns[col])
429429+ for (row = 0; row<nOfRows; row++)
430430+ distMatrix[row + nOfRows*col] -= h;
431431+432432+ /* move to step 3 */
433433+ step3(assignment, distMatrix, starMatrix, newStarMatrix, primeMatrix, coveredColumns, coveredRows, nOfRows, nOfColumns, minDim);
434434+}
+23
src/external/hungarian/LICENSE
···11+Copyright (c) 2016, mcximing
22+All rights reserved.
33+44+Redistribution and use in source and binary forms, with or without
55+modification, are permitted provided that the following conditions are met:
66+77+* Redistributions of source code must retain the above copyright notice, this
88+ list of conditions and the following disclaimer.
99+1010+* Redistributions in binary form must reproduce the above copyright notice,
1111+ this list of conditions and the following disclaimer in the documentation
1212+ and/or other materials provided with the distribution.
1313+1414+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
1515+AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
1616+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
1717+DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
1818+FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
1919+DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
2020+SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
2121+CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
2222+OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
2323+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+7
src/external/hungarian/README.txt
···11+22+Hungarian algorithm, also known as Munkres algorithm or Kuhn-Munkres algorithm, is a method for solving the assignment problem, for example assigning workers to jobs, which goal is to compute the optimal assignment that minimizes the total cost, and the like.
33+44+This is a C++ wrapper with slight modification of a hungarian algorithm implementation by Markus Buehren.
55+The original code is a few mex-functions for use in MATLAB, found here:
66+http://www.mathworks.com/matlabcentral/fileexchange/6543-functions-for-the-rectangular-assignment-problem
77+