多目标跟踪之数据关联算法——匈牙利算法

零、Track和Detection的cost matrix，distance metric。距离计算的方式有如下几种：

1.  距离cost distance metric，track和detection的距离矩阵。
   外观距离appearance distance，来自检测切片ROI的网络特征提取；——余弦距离
   运动模型距离 马氏距离，来自检测-跟踪的kalman校正距离，马氏距离==二维高斯分布，利用相关可以去除维度间的量纲影响，并考虑维度间的相关性。比如x,y,a,h四者的相关性和量纲差异，可以不被考虑。
2. 距离矩阵的分配算法——匈牙利算法Hungarian，这里指的是Kuhn&Munkres Hungarian， KM匈牙利算法（加权的匈牙利算法）
   ① KM匈牙利，n个轨迹，m个检测。padding方式为max(n,m)
   ② KM匈牙利，直接对矩形n×m进行分配，不进行任何形式的padding
   ③ KM匈牙利，padding方式为(n+m)*(n+m)，主对角有数值。其余分别为inf和0
   这些算法的优缺点和差异分别为什么？
3. 距离矩阵进入匈牙利算法可以有不同的padding形式，他们对分配结果的差异是怎样的？
   
   

track_t CTrack::CalcDistCenter(const CRegion& reg) const
{
    Point_t diff = m_predictionPoint - reg.m_rrect.center;
    return sqrtf(sqr(diff.x) + sqr(diff.y));
}

track_t CTrack::CalcDistRect(const CRegion& reg) const
{
    std::array<track_t, 5> diff;
    diff[0] = reg.m_rrect.center.x - m_lastRegion.m_rrect.center.x;
    diff[1] = reg.m_rrect.center.y - m_lastRegion.m_rrect.center.y;
    diff[2] = static_cast<track_t>(m_lastRegion.m_rrect.size.width - reg.m_rrect.size.width);
    diff[3] = static_cast<track_t>(m_lastRegion.m_rrect.size.height - reg.m_rrect.size.height);
    diff[4] = static_cast<track_t>(m_lastRegion.m_rrect.angle - reg.m_rrect.angle);

    track_t dist = 0;
    for (size_t i = 0; i < diff.size(); ++i)
    {
        dist += sqr(diff[i]);
    }
    return sqrtf(dist);
}

track_t CTrack::CalcDistJaccard(const CRegion& reg) const
{
    track_t intArea = static_cast<track_t>((reg.m_brect & m_lastRegion.m_brect).area());
    track_t unionArea = static_cast<track_t>(reg.m_brect.area() + m_lastRegion.m_brect.area() - intArea) + 1e-6;

    return std::fabs(1 - intArea / unionArea);
}
track_t CTrack::CalcDistHist(const RegionEmbedding& embedding) const
{
	track_t res = 1;

    if (!embedding.m_hist.empty() && !m_regionEmbedding.m_hist.empty())
	{
#if (((CV_VERSION_MAJOR == 4) && (CV_VERSION_MINOR < 1)) || (CV_VERSION_MAJOR == 3))
		res = static_cast<track_t>(cv::compareHist(embedding.m_hist, m_regionEmbedding.m_hist, CV_COMP_BHATTACHARYYA));
        //res = 1.f - static_cast<track_t>(cv::compareHist(hist, m_regionEmbedding.m_hist, CV_COMP_CORREL));
#else
        res = static_cast<track_t>(cv::compareHist(embedding.m_hist, m_regionEmbedding.m_hist, cv::HISTCMP_BHATTACHARYYA));
#endif
	}
    else
    {
        assert(0);
        CV_Assert(!embedding.m_hist.empty());
        CV_Assert(!m_regionEmbedding.m_hist.empty());
    }
	return res;
}
track_t CTrack::CalcMahalanobisDist(const cv::RotatedRect& rrect) const
{
    cv::Mat res1, predictPoint;
    // res1 = Hn * Pn+1|n+1 * Hn^T + Rn+1	error covariance
    // res2 = Hn * Xn+1|n
    m_kalman.GetPtStateAndResCov(res1, predictPoint);

    double mahaDist = 0.0;
    if (!res1.empty() && !predictPoint.empty())
    {
        cv::Mat icovar_Pn;
        cv::invert(res1, icovar_Pn, cv::DECOMP_SVD);
        cv::Mat measurePoint;
        if (predictPoint.rows == 2)			// PointUpdate
            measurePoint = (cv::Mat_<track_t>(2, 1) << rrect.center.x, rrect.center.y); // detection
        else
            measurePoint = (cv::Mat_<track_t>(4, 1) << rrect.center.x, rrect.center.y, rrect.size.width, rrect.size.height); // predict
        mahaDist = cv::Mahalanobis(measurePoint, predictPoint, icovar_Pn);
        mahaDist += std::log(cv::determinant(res1));
    }
    return static_cast<track_t>(mahaDist);
}
std::pair<track_t, bool> CTrack::CalcCosine(const RegionEmbedding& embedding) const
{
	track_t res = 1;
	if (!embedding.m_embedding.empty() && !m_regionEmbedding.m_embedding.empty())
	{
		double xy = embedding.m_embedding.dot(m_regionEmbedding.m_embedding);
		double norm = sqrt(embedding.m_embDot * m_regionEmbedding.m_embDot) + 1e-6;
#if 0
        res = 1.f - 0.5f * fabs(static_cast<float>(xy / norm));
#else
        res = 0.5f * static_cast<float>(1.0 - xy / norm);
#endif
        //std::cout << "CTrack::CalcCosine: " << embedding.m_embedding.size() << " - " << m_regionEmbedding.m_embedding.size() << " = " << res << std::endl;
        return { res, true };
	}
    else
    {
        //assert(0);
        //CV_Assert(!embedding.m_embedding.empty());
        //CV_Assert(!m_regionEmbedding.m_embedding.empty());
        return { 0, false };
    }
}

---

 一、匈牙利算法维基百科：https://zh.wikipedia.org/wiki/匈牙利算法

• 匈牙利算法最早是由匈牙利数学家Dénes Kőnig康尼格和Jenő Egerváry用来求矩阵中0元素的个数的一种方法。由此他证明『矩阵中独立0元素的最多个数等于能覆盖cover所有0元素的最少直线数』
• 匈牙利算法是一种在多项式时间内求解任务分配问题的组合优化算法，并推动了后来的原始对偶方法。
• 1955年，美国数学家哈罗德·库恩W.W.Kuhn（库恩）Harold W. Kuhn在求解著名的指派问题时引用了这一结论，并对具体算法做了改进，仍然称为匈牙利算法。此算法之所以被称作匈牙利算法，是因为算法很大一部分是基于以前匈牙利数学家Dénes Kőnig和Jenő Egerváry的工作之上创建起来的。[1][2]
• 1957年，詹姆士·芒克勒斯Munkres在1957年回顾了该算法，并发现它的时间复杂度为（强）多项式时间。[3] 此后该算法被称为Kuhn–Munkres算法或Munkres分配算法。原始算法的时间复杂度为 ，但杰克·爱德蒙斯与卡普发现可以修改算法达到 运行时间，富泽也独立发现了这一点。L·R·福特和D·R·福尔克森将该方法推广到了一般运输问题。2006年发现卡尔·雅可比 Jacobis bound在19世纪就解决了指派问题，该解法在他死后在1890年以拉丁文发表。[4]

参考书目

• R.E. Burkard, M. Dell'Amico, S. Martello: Assignment Problems (Revised reprint). SIAM, Philadelphia (PA.) 2012. ISBN 978-1-61197-222-1
• M. Fischetti, "Lezioni di Ricerca Operativa", Edizioni Libreria Progetto Padova, Italia, 1995.
• R. Ahuja, T. Magnanti, J. Orlin, "Network Flows", Prentice Hall, 1993.
• S. Martello, "Jeno Egerváry: from the origins of the Hungarian algorithm to satellite communication". Central European Journal of Operations Research 18, 47–58, 2010

参考文献

1. ^ Harold W. Kuhn, "The Hungarian Method for the assignment problem", Naval Research Logistics Quarterly, 2: 83–97, 1955. Kuhn's original publication.
2. ^ Harold W. Kuhn, "Variants of the Hungarian method for assignment problems", Naval Research Logistics Quarterly, 3: 253–258, 1956.
    https://econweb.ucsd.edu/~v2crawford/hungar.pdf  Kuhn的paper1960年9月29日。
3. ^ J. Munkres, "Algorithms for the Assignment and Transportation Problems", Journal of the Society for Industrial and Applied Mathematics, 5(1):32–38, 1957 March.
4. ^ JACOBI'S BOUND

外部链接

• Bruff, Derek, "The Assignment Problem and the Hungarian Method", [1] （页面存档备份，存于互联网档案馆）
• Mordecai J. Golin, Bipartite Matching and the Hungarian Method, Course Notes, Hong Kong University of Science and Technology.
• R. A. Pilgrim, Munkres' Assignment Algorithm. Modified for Rectangular Matrices （页面存档备份，存于互联网档案馆）, Course notes, Murray State University.
• Mike Dawes, The Optimal Assignment Problem, Course notes, University of Western Ontario.
• On Kuhn's Hungarian Method – A tribute from Hungary, András Frank, Egervary Research Group, Pazmany P. setany 1/C, H1117, Budapest, Hungary.
• Lecture: Fundamentals of Operations Research - Assignment Problem - Hungarian Algorithm, Prof. G. Srinivasan, Department of Management Studies, IIT Madras.
• Extension: Assignment sensitivity analysis (with O(n^4) time complexity), Liu, Shell.
• Solve any Assignment Problem online, provides a step by step explanation of the Hungarian Algorithm.

各种语言的算法实现链接

（请注意，并非所有这些都满足 时间约束。） 

• C implementation with O ( n 3 )  time complexity
• Java implementation of O ( n 3 ) time variant
• Python implementation (see also here)
• Ruby implementation with unit tests
• C# implementation
• D implementation with unit tests (port of the Java O ( n 3 )  version)
• Online interactive implementation Please note that this implements a variant of the algorithm as described above.
• Graphical implementation with options (Java applet)
• Serial and parallel implementations.
• Implementation in Matlab and C
• Perl implementation
• Lisp implementation