/*M/////////////////////////////////////////////////////////////////////////////////////// // // IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING. // // By downloading, copying, installing or using the software you agree to this license. // If you do not agree to this license, do not download, install, // copy or use the software. // // // License Agreement // For Open Source Computer Vision Library // // Copyright (C) 2000-2008, Intel Corporation, all rights reserved. // Copyright (C) 2009, Willow Garage Inc., all rights reserved. // Third party copyrights are property of their respective owners. // // Redistribution and use in source and binary forms, with or without modification, // are permitted provided that the following conditions are met: // // * Redistribution's of source code must retain the above copyright notice, // this list of conditions and the following disclaimer. // // * Redistribution's in binary form must reproduce the above copyright notice, // this list of conditions and the following disclaimer in the documentation // and/or other materials provided with the distribution. // // * The name of the copyright holders may not be used to endorse or promote products // derived from this software without specific prior written permission. // // This software is provided by the copyright holders and contributors "as is" and // any express or implied warranties, including, but not limited to, the implied // warranties of merchantability and fitness for a particular purpose are disclaimed. // In no event shall the Intel Corporation or contributors be liable for any direct, // indirect, incidental, special, exemplary, or consequential damages // (including, but not limited to, procurement of substitute goods or services; // loss of use, data, or profits; or business interruption) however caused // and on any theory of liability, whether in contract, strict liability, // or tort (including negligence or otherwise) arising in any way out of // the use of this software, even if advised of the possibility of such damage. // //M*/ /* * Implementation of an optimized EMD for histograms based in * the papers "EMD-L1: An efficient and Robust Algorithm * for comparing histogram-based descriptors", by Haibin Ling and * Kazunori Okuda; and "The Earth Mover's Distance is the Mallows * Distance: Some Insights from Statistics", by Elizaveta Levina and * Peter Bickel, based on HAIBIN LING AND KAZUNORI OKADA implementation. */ #include "precomp.hpp" #include "emdL1_def.hpp" #include <limits> /****************************************************************************************\ * EMDL1 Class * \****************************************************************************************/ float EmdL1::getEMDL1(cv::Mat &sig1, cv::Mat &sig2) { // Initialization CV_Assert((sig1.rows==sig2.rows) && (sig1.cols==sig2.cols) && (!sig1.empty()) && (!sig2.empty())); if(!initBaseTrees(sig1.rows, 1)) return -1; float *H1=new float[sig1.rows], *H2 = new float[sig2.rows]; for (int ii=0; ii<sig1.rows; ii++) { H1[ii]=sig1.at<float>(ii,0); H2[ii]=sig2.at<float>(ii,0); } fillBaseTrees(H1,H2); // Initialize histograms greedySolution(); // Construct an initial Basic Feasible solution initBVTree(); // Initialize BVTree // Iteration bool bOptimal = false; m_nItr = 0; while(!bOptimal && m_nItr<nMaxIt) { // Derive U=(u_ij) for row i and column j if(m_nItr==0) updateSubtree(m_pRoot); else updateSubtree(m_pEnter->pChild); // Optimality test bOptimal = isOptimal(); // Find new solution if(!bOptimal) findNewSolution(); ++m_nItr; } delete [] H1; delete [] H2; // Output the total flow return compuTotalFlow(); } void EmdL1::setMaxIteration(int _nMaxIt) { nMaxIt=_nMaxIt; } //-- SubFunctions called in the EMD algorithm bool EmdL1::initBaseTrees(int n1, int n2, int n3) { if(binsDim1==n1 && binsDim2==n2 && binsDim3==n3) return true; binsDim1 = n1; binsDim2 = n2; binsDim3 = n3; if(binsDim1==0 || binsDim2==0) dimension = 0; else dimension = (binsDim3==0)?2:3; if(dimension==2) { m_Nodes.resize(binsDim1); m_EdgesUp.resize(binsDim1); m_EdgesRight.resize(binsDim1); for(int i1=0; i1<binsDim1; i1++) { m_Nodes[i1].resize(binsDim2); m_EdgesUp[i1].resize(binsDim2); m_EdgesRight[i1].resize(binsDim2); } m_NBVEdges.resize(binsDim1*binsDim2*4+2); m_auxQueue.resize(binsDim1*binsDim2+2); m_fromLoop.resize(binsDim1*binsDim2+2); m_toLoop.resize(binsDim1*binsDim2+2); } else if(dimension==3) { m_3dNodes.resize(binsDim1); m_3dEdgesUp.resize(binsDim1); m_3dEdgesRight.resize(binsDim1); m_3dEdgesDeep.resize(binsDim1); for(int i1=0; i1<binsDim1; i1++) { m_3dNodes[i1].resize(binsDim2); m_3dEdgesUp[i1].resize(binsDim2); m_3dEdgesRight[i1].resize(binsDim2); m_3dEdgesDeep[i1].resize(binsDim2); for(int i2=0; i2<binsDim2; i2++) { m_3dNodes[i1][i2].resize(binsDim3); m_3dEdgesUp[i1][i2].resize(binsDim3); m_3dEdgesRight[i1][i2].resize(binsDim3); m_3dEdgesDeep[i1][i2].resize(binsDim3); } } m_NBVEdges.resize(binsDim1*binsDim2*binsDim3*6+4); m_auxQueue.resize(binsDim1*binsDim2*binsDim3+4); m_fromLoop.resize(binsDim1*binsDim2*binsDim3+4); m_toLoop.resize(binsDim1*binsDim2*binsDim3+2); } else return false; return true; } bool EmdL1::fillBaseTrees(float *H1, float *H2) { //- Set global counters m_pRoot = NULL; // Graph initialization float *p1 = H1; float *p2 = H2; if(dimension==2) { for(int c=0; c<binsDim2; c++) { for(int r=0; r<binsDim1; r++) { //- initialize nodes and links m_Nodes[r][c].pos[0] = r; m_Nodes[r][c].pos[1] = c; m_Nodes[r][c].d = *(p1++)-*(p2++); m_Nodes[r][c].pParent = NULL; m_Nodes[r][c].pChild = NULL; m_Nodes[r][c].iLevel = -1; //- initialize edges // to the right m_EdgesRight[r][c].pParent = &(m_Nodes[r][c]); m_EdgesRight[r][c].pChild = &(m_Nodes[r][(c+1)%binsDim2]); m_EdgesRight[r][c].flow = 0; m_EdgesRight[r][c].iDir = 1; m_EdgesRight[r][c].pNxt = NULL; // to the upward m_EdgesUp[r][c].pParent = &(m_Nodes[r][c]); m_EdgesUp[r][c].pChild = &(m_Nodes[(r+1)%binsDim1][c]); m_EdgesUp[r][c].flow = 0; m_EdgesUp[r][c].iDir = 1; m_EdgesUp[r][c].pNxt = NULL; } } } else if(dimension==3) { for(int z=0; z<binsDim3; z++) { for(int c=0; c<binsDim2; c++) { for(int r=0; r<binsDim1; r++) { //- initialize nodes and edges m_3dNodes[r][c][z].pos[0] = r; m_3dNodes[r][c][z].pos[1] = c; m_3dNodes[r][c][z].pos[2] = z; m_3dNodes[r][c][z].d = *(p1++)-*(p2++); m_3dNodes[r][c][z].pParent = NULL; m_3dNodes[r][c][z].pChild = NULL; m_3dNodes[r][c][z].iLevel = -1; //- initialize edges // to the upward m_3dEdgesUp[r][c][z].pParent= &(m_3dNodes[r][c][z]); m_3dEdgesUp[r][c][z].pChild = &(m_3dNodes[(r+1)%binsDim1][c][z]); m_3dEdgesUp[r][c][z].flow = 0; m_3dEdgesUp[r][c][z].iDir = 1; m_3dEdgesUp[r][c][z].pNxt = NULL; // to the right m_3dEdgesRight[r][c][z].pParent = &(m_3dNodes[r][c][z]); m_3dEdgesRight[r][c][z].pChild = &(m_3dNodes[r][(c+1)%binsDim2][z]); m_3dEdgesRight[r][c][z].flow = 0; m_3dEdgesRight[r][c][z].iDir = 1; m_3dEdgesRight[r][c][z].pNxt = NULL; // to the deep m_3dEdgesDeep[r][c][z].pParent = &(m_3dNodes[r][c][z]); m_3dEdgesDeep[r][c][z].pChild = &(m_3dNodes[r][c])[(z+1)%binsDim3]; m_3dEdgesDeep[r][c][z].flow = 0; m_3dEdgesDeep[r][c][z].iDir = 1; m_3dEdgesDeep[r][c][z].pNxt = NULL; } } } } return true; } bool EmdL1::greedySolution() { return dimension==2?greedySolution2():greedySolution3(); } bool EmdL1::greedySolution2() { //- Prepare auxiliary array, D=H1-H2 int c,r; floatArray2D D(binsDim1); for(r=0; r<binsDim1; r++) { D[r].resize(binsDim2); for(c=0; c<binsDim2; c++) D[r][c] = m_Nodes[r][c].d; } // compute integrated values along each dimension std::vector<float> d2s(binsDim2); d2s[0] = 0; for(c=0; c<binsDim2-1; c++) { d2s[c+1] = d2s[c]; for(r=0; r<binsDim1; r++) d2s[c+1]-= D[r][c]; } std::vector<float> d1s(binsDim1); d1s[0] = 0; for(r=0; r<binsDim1-1; r++) { d1s[r+1] = d1s[r]; for(c=0; c<binsDim2; c++) d1s[r+1]-= D[r][c]; } //- Greedy algorithm for initial solution cvPEmdEdge pBV; float dFlow; bool bUpward = false; nNBV = 0; // number of NON-BV edges for(c=0; c<binsDim2-1; c++) for(r=0; r<binsDim1; r++) { dFlow = D[r][c]; bUpward = (r<binsDim1-1) && (fabs(dFlow+d2s[c+1]) > fabs(dFlow+d1s[r+1])); // Move upward or right // modify basic variables, record BV and related values if(bUpward) { // move to up pBV = &(m_EdgesUp[r][c]); m_NBVEdges[nNBV++] = &(m_EdgesRight[r][c]); D[r+1][c] += dFlow; // auxiliary matrix maintenance d1s[r+1] += dFlow; // auxiliary matrix maintenance } else { // move to right, no other choice pBV = &(m_EdgesRight[r][c]); if(r<binsDim1-1) m_NBVEdges[nNBV++] = &(m_EdgesUp[r][c]); D[r][c+1] += dFlow; // auxiliary matrix maintenance d2s[c+1] += dFlow; // auxiliary matrix maintenance } pBV->pParent->pChild = pBV; pBV->flow = fabs(dFlow); pBV->iDir = dFlow>0; // 1:outward, 0:inward } //- rightmost column, no choice but move upward c = binsDim2-1; for(r=0; r<binsDim1-1; r++) { dFlow = D[r][c]; pBV = &(m_EdgesUp[r][c]); D[r+1][c] += dFlow; // auxiliary matrix maintenance pBV->pParent->pChild= pBV; pBV->flow = fabs(dFlow); pBV->iDir = dFlow>0; // 1:outward, 0:inward } return true; } bool EmdL1::greedySolution3() { //- Prepare auxiliary array, D=H1-H2 int i1,i2,i3; std::vector<floatArray2D> D(binsDim1); for(i1=0; i1<binsDim1; i1++) { D[i1].resize(binsDim2); for(i2=0; i2<binsDim2; i2++) { D[i1][i2].resize(binsDim3); for(i3=0; i3<binsDim3; i3++) D[i1][i2][i3] = m_3dNodes[i1][i2][i3].d; } } // compute integrated values along each dimension std::vector<float> d1s(binsDim1); d1s[0] = 0; for(i1=0; i1<binsDim1-1; i1++) { d1s[i1+1] = d1s[i1]; for(i2=0; i2<binsDim2; i2++) { for(i3=0; i3<binsDim3; i3++) d1s[i1+1] -= D[i1][i2][i3]; } } std::vector<float> d2s(binsDim2); d2s[0] = 0; for(i2=0; i2<binsDim2-1; i2++) { d2s[i2+1] = d2s[i2]; for(i1=0; i1<binsDim1; i1++) { for(i3=0; i3<binsDim3; i3++) d2s[i2+1] -= D[i1][i2][i3]; } } std::vector<float> d3s(binsDim3); d3s[0] = 0; for(i3=0; i3<binsDim3-1; i3++) { d3s[i3+1] = d3s[i3]; for(i1=0; i1<binsDim1; i1++) { for(i2=0; i2<binsDim2; i2++) d3s[i3+1] -= D[i1][i2][i3]; } } //- Greedy algorithm for initial solution cvPEmdEdge pBV; float dFlow, f1,f2,f3; nNBV = 0; // number of NON-BV edges for(i3=0; i3<binsDim3; i3++) { for(i2=0; i2<binsDim2; i2++) { for(i1=0; i1<binsDim1; i1++) { if(i3==binsDim3-1 && i2==binsDim2-1 && i1==binsDim1-1) break; //- determine which direction to move, either right or upward dFlow = D[i1][i2][i3]; f1 = (i1<(binsDim1-1))?fabs(dFlow+d1s[i1+1]):std::numeric_limits<float>::max(); f2 = (i2<(binsDim2-1))?fabs(dFlow+d2s[i2+1]):std::numeric_limits<float>::max(); f3 = (i3<(binsDim3-1))?fabs(dFlow+d3s[i3+1]):std::numeric_limits<float>::max(); if(f1<f2 && f1<f3) { pBV = &(m_3dEdgesUp[i1][i2][i3]); // up if(i2<binsDim2-1) m_NBVEdges[nNBV++] = &(m_3dEdgesRight[i1][i2][i3]); // right if(i3<binsDim3-1) m_NBVEdges[nNBV++] = &(m_3dEdgesDeep[i1][i2][i3]); // deep D[i1+1][i2][i3] += dFlow; // maintain auxiliary matrix d1s[i1+1] += dFlow; } else if(f2<f3) { pBV = &(m_3dEdgesRight[i1][i2][i3]); // right if(i1<binsDim1-1) m_NBVEdges[nNBV++] = &(m_3dEdgesUp[i1][i2][i3]); // up if(i3<binsDim3-1) m_NBVEdges[nNBV++] = &(m_3dEdgesDeep[i1][i2][i3]); // deep D[i1][i2+1][i3] += dFlow; // maintain auxiliary matrix d2s[i2+1] += dFlow; } else { pBV = &(m_3dEdgesDeep[i1][i2][i3]); // deep if(i2<binsDim2-1) m_NBVEdges[nNBV++] = &(m_3dEdgesRight[i1][i2][i3]); // right if(i1<binsDim1-1) m_NBVEdges[nNBV++] = &(m_3dEdgesUp[i1][i2][i3]); // up D[i1][i2][i3+1] += dFlow; // maintain auxiliary matrix d3s[i3+1] += dFlow; } pBV->flow = fabs(dFlow); pBV->iDir = dFlow>0; // 1:outward, 0:inward pBV->pParent->pChild= pBV; } } } return true; } void EmdL1::initBVTree() { // initialize BVTree from the initial BF solution //- Using the center of the graph as the root int r = (int)(0.5*binsDim1-.5); int c = (int)(0.5*binsDim2-.5); int z = (int)(0.5*binsDim3-.5); m_pRoot = dimension==2 ? &(m_Nodes[r][c]) : &(m_3dNodes[r][c][z]); m_pRoot->u = 0; m_pRoot->iLevel = 0; m_pRoot->pParent= NULL; m_pRoot->pPEdge = NULL; //- Prepare a queue m_auxQueue[0] = m_pRoot; int nQueue = 1; // length of queue int iQHead = 0; // head of queue //- Recursively build subtrees cvPEmdEdge pCurE=NULL, pNxtE=NULL; cvPEmdNode pCurN=NULL, pNxtN=NULL; int nBin = binsDim1*binsDim2*std::max(binsDim3,1); while(iQHead<nQueue && nQueue<nBin) { pCurN = m_auxQueue[iQHead++]; // pop out from queue r = pCurN->pos[0]; c = pCurN->pos[1]; z = pCurN->pos[2]; // check connection from itself pCurE = pCurN->pChild; // the initial child from initial solution if(pCurE) { pNxtN = pCurE->pChild; pNxtN->pParent = pCurN; pNxtN->pPEdge = pCurE; m_auxQueue[nQueue++] = pNxtN; } // check four neighbor nodes int nNB = dimension==2?4:6; for(int k=0;k<nNB;k++) { if(dimension==2) { if(k==0 && c>0) pNxtN = &(m_Nodes[r][c-1]); // left else if(k==1 && r>0) pNxtN = &(m_Nodes[r-1][c]); // down else if(k==2 && c<binsDim2-1) pNxtN = &(m_Nodes[r][c+1]); // right else if(k==3 && r<binsDim1-1) pNxtN = &(m_Nodes[r+1][c]); // up else continue; } else if(dimension==3) { if(k==0 && c>0) pNxtN = &(m_3dNodes[r][c-1][z]); // left else if(k==1 && c<binsDim2-1) pNxtN = &(m_3dNodes[r][c+1][z]); // right else if(k==2 && r>0) pNxtN = &(m_3dNodes[r-1][c][z]); // down else if(k==3 && r<binsDim1-1) pNxtN = &(m_3dNodes[r+1][c][z]); // up else if(k==4 && z>0) pNxtN = &(m_3dNodes[r][c][z-1]); // shallow else if(k==5 && z<binsDim3-1) pNxtN = &(m_3dNodes[r][c][z+1]); // deep else continue; } if(pNxtN != pCurN->pParent) { CV_Assert(pNxtN != NULL); pNxtE = pNxtN->pChild; if(pNxtE && pNxtE->pChild==pCurN) // has connection { pNxtN->pParent = pCurN; pNxtN->pPEdge = pNxtE; pNxtN->pChild = NULL; m_auxQueue[nQueue++] = pNxtN; pNxtE->pParent = pCurN; // reverse direction pNxtE->pChild = pNxtN; pNxtE->iDir = !pNxtE->iDir; if(pCurE) pCurE->pNxt = pNxtE; // add to edge list else pCurN->pChild = pNxtE; pCurE = pNxtE; } } } } } void EmdL1::updateSubtree(cvPEmdNode pRoot) { // Initialize auxiliary queue m_auxQueue[0] = pRoot; int nQueue = 1; // queue length int iQHead = 0; // head of queue // BFS browing cvPEmdNode pCurN=NULL,pNxtN=NULL; cvPEmdEdge pCurE=NULL; while(iQHead<nQueue) { pCurN = m_auxQueue[iQHead++]; // pop out from queue pCurE = pCurN->pChild; // browsing all children while(pCurE) { pNxtN = pCurE->pChild; pNxtN->iLevel = pCurN->iLevel+1; pNxtN->u = pCurE->iDir ? (pCurN->u - 1) : (pCurN->u + 1); pCurE = pCurE->pNxt; m_auxQueue[nQueue++] = pNxtN; } } } bool EmdL1::isOptimal() { int iC, iMinC = 0; cvPEmdEdge pE; m_pEnter = NULL; m_iEnter = -1; // test each NON-BV edges for(int k=0; k<nNBV; ++k) { pE = m_NBVEdges[k]; iC = 1 - pE->pParent->u + pE->pChild->u; if(iC<iMinC) { iMinC = iC; m_iEnter= k; } else { // Try reversing the direction iC = 1 + pE->pParent->u - pE->pChild->u; if(iC<iMinC) { iMinC = iC; m_iEnter= k; } } } if(m_iEnter>=0) { m_pEnter = m_NBVEdges[m_iEnter]; if(iMinC == (1 - m_pEnter->pChild->u + m_pEnter->pParent->u)) { // reverse direction cvPEmdNode pN = m_pEnter->pParent; m_pEnter->pParent = m_pEnter->pChild; m_pEnter->pChild = pN; } m_pEnter->iDir = 1; } return m_iEnter==-1; } void EmdL1::findNewSolution() { // Find loop formed by adding the Enter BV edge. findLoopFromEnterBV(); // Modify flow values along the loop cvPEmdEdge pE = NULL; CV_Assert(m_pLeave != NULL); float minFlow = m_pLeave->flow; int k; for(k=0; k<m_iFrom; k++) { pE = m_fromLoop[k]; if(pE->iDir) pE->flow += minFlow; // outward else pE->flow -= minFlow; // inward } for(k=0; k<m_iTo; k++) { pE = m_toLoop[k]; if(pE->iDir) pE->flow -= minFlow; // outward else pE->flow += minFlow; // inward } // Update BV Tree, removing the Leaving-BV edge cvPEmdNode pLParentN = m_pLeave->pParent; cvPEmdNode pLChildN = m_pLeave->pChild; cvPEmdEdge pPreE = pLParentN->pChild; if(pPreE==m_pLeave) { pLParentN->pChild = m_pLeave->pNxt; // Leaving-BV is the first child } else { while(pPreE->pNxt != m_pLeave) pPreE = pPreE->pNxt; pPreE->pNxt = m_pLeave->pNxt; // remove Leaving-BV from child list } pLChildN->pParent = NULL; pLChildN->pPEdge = NULL; m_NBVEdges[m_iEnter]= m_pLeave; // put the leaving-BV into the NBV array // Add the Enter BV edge cvPEmdNode pEParentN = m_pEnter->pParent; cvPEmdNode pEChildN = m_pEnter->pChild; m_pEnter->flow = minFlow; m_pEnter->pNxt = pEParentN->pChild; // insert the Enter BV as the first child pEParentN->pChild = m_pEnter; // of its parent // Recursively update the tree start from pEChildN cvPEmdNode pPreN = pEParentN; cvPEmdNode pCurN = pEChildN; cvPEmdNode pNxtN; cvPEmdEdge pNxtE, pPreE0; pPreE = m_pEnter; while(pCurN) { pNxtN = pCurN->pParent; pNxtE = pCurN->pPEdge; pCurN->pParent = pPreN; pCurN->pPEdge = pPreE; if(pNxtN) { // remove the edge from pNxtN's child list if(pNxtN->pChild==pNxtE) { pNxtN->pChild = pNxtE->pNxt; // first child } else { pPreE0 = pNxtN->pChild; while(pPreE0->pNxt != pNxtE) pPreE0 = pPreE0->pNxt; pPreE0->pNxt = pNxtE->pNxt; // remove Leaving-BV from child list } // reverse the parent-child direction pNxtE->pParent = pCurN; pNxtE->pChild = pNxtN; pNxtE->iDir = !pNxtE->iDir; pNxtE->pNxt = pCurN->pChild; pCurN->pChild = pNxtE; pPreE = pNxtE; pPreN = pCurN; } pCurN = pNxtN; } // Update U at the child of the Enter BV pEChildN->u = m_pEnter->iDir?(pEParentN->u-1):(pEParentN->u + 1); pEChildN->iLevel = pEParentN->iLevel+1; } void EmdL1::findLoopFromEnterBV() { // Initialize Leaving-BV edge float minFlow = std::numeric_limits<float>::max(); cvPEmdEdge pE = NULL; int iLFlag = 0; // 0: in the FROM list, 1: in the TO list // Using two loop list to store the loop nodes cvPEmdNode pFrom = m_pEnter->pParent; cvPEmdNode pTo = m_pEnter->pChild; m_iFrom = 0; m_iTo = 0; m_pLeave = NULL; // Trace back to make pFrom and pTo at the same level while(pFrom->iLevel > pTo->iLevel) { pE = pFrom->pPEdge; m_fromLoop[m_iFrom++] = pE; if(!pE->iDir && pE->flow<minFlow) { minFlow = pE->flow; m_pLeave = pE; iLFlag = 0; // 0: in the FROM list } pFrom = pFrom->pParent; } while(pTo->iLevel > pFrom->iLevel) { pE = pTo->pPEdge; m_toLoop[m_iTo++] = pE; if(pE->iDir && pE->flow<minFlow) { minFlow = pE->flow; m_pLeave = pE; iLFlag = 1; // 1: in the TO list } pTo = pTo->pParent; } // Trace pTo and pFrom simultaneously till find their common ancester while(pTo!=pFrom) { pE = pFrom->pPEdge; m_fromLoop[m_iFrom++] = pE; if(!pE->iDir && pE->flow<minFlow) { minFlow = pE->flow; m_pLeave = pE; iLFlag = 0; // 0: in the FROM list, 1: in the TO list } pFrom = pFrom->pParent; pE = pTo->pPEdge; m_toLoop[m_iTo++] = pE; if(pE->iDir && pE->flow<minFlow) { minFlow = pE->flow; m_pLeave = pE; iLFlag = 1; // 0: in the FROM list, 1: in the TO list } pTo = pTo->pParent; } // Reverse the direction of the Enter BV edge if necessary if(iLFlag==0) { cvPEmdNode pN = m_pEnter->pParent; m_pEnter->pParent = m_pEnter->pChild; m_pEnter->pChild = pN; m_pEnter->iDir = !m_pEnter->iDir; } } float EmdL1::compuTotalFlow() { // Computing the total flow as the final distance float f = 0; // Initialize auxiliary queue m_auxQueue[0] = m_pRoot; int nQueue = 1; // length of queue int iQHead = 0; // head of queue // BFS browing the tree cvPEmdNode pCurN=NULL,pNxtN=NULL; cvPEmdEdge pCurE=NULL; while(iQHead<nQueue) { pCurN = m_auxQueue[iQHead++]; // pop out from queue pCurE = pCurN->pChild; // browsing all children while(pCurE) { f += pCurE->flow; pNxtN = pCurE->pChild; pCurE = pCurE->pNxt; m_auxQueue[nQueue++] = pNxtN; } } return f; } /****************************************************************************************\ * EMDL1 Function * \****************************************************************************************/ float cv::EMDL1(InputArray _signature1, InputArray _signature2) { CV_INSTRUMENT_REGION(); Mat signature1 = _signature1.getMat(), signature2 = _signature2.getMat(); EmdL1 emdl1; return emdl1.getEMDL1(signature1, signature2); }