1 /*M/////////////////////////////////////////////////////////////////////////////////////// 2 // 3 // IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING. 4 // 5 // By downloading, copying, installing or using the software you agree to this license. 6 // If you do not agree to this license, do not download, install, 7 // copy or use the software. 8 // 9 // 10 // License Agreement 11 // For Open Source Computer Vision Library 12 // 13 // Copyright (C) 2000-2008, Intel Corporation, all rights reserved. 14 // Copyright (C) 2009, Willow Garage Inc., all rights reserved. 15 // Third party copyrights are property of their respective owners. 16 // 17 // Redistribution and use in source and binary forms, with or without modification, 18 // are permitted provided that the following conditions are met: 19 // 20 // * Redistribution's of source code must retain the above copyright notice, 21 // this list of conditions and the following disclaimer. 22 // 23 // * Redistribution's in binary form must reproduce the above copyright notice, 24 // this list of conditions and the following disclaimer in the documentation 25 // and/or other materials provided with the distribution. 26 // 27 // * The name of the copyright holders may not be used to endorse or promote products 28 // derived from this software without specific prior written permission. 29 // 30 // This software is provided by the copyright holders and contributors "as is" and 31 // any express or implied warranties, including, but not limited to, the implied 32 // warranties of merchantability and fitness for a particular purpose are disclaimed. 33 // In no event shall the Intel Corporation or contributors be liable for any direct, 34 // indirect, incidental, special, exemplary, or consequential damages 35 // (including, but not limited to, procurement of substitute goods or services; 36 // loss of use, data, or profits; or business interruption) however caused 37 // and on any theory of liability, whether in contract, strict liability, 38 // or tort (including negligence or otherwise) arising in any way out of 39 // the use of this software, even if advised of the possibility of such damage. 40 // 41 //M*/ 42 43 /* 44 * Implementation of an optimized EMD for histograms based in 45 * the papers "EMD-L1: An efficient and Robust Algorithm 46 * for comparing histogram-based descriptors", by Haibin Ling and 47 * Kazunori Okuda; and "The Earth Mover's Distance is the Mallows 48 * Distance: Some Insights from Statistics", by Elizaveta Levina and 49 * Peter Bickel, based on HAIBIN LING AND KAZUNORI OKADA implementation. 50 */ 51 52 #include "precomp.hpp" 53 #include "emdL1_def.hpp" 54 #include <limits> 55 56 /****************************************************************************************\ 57 * EMDL1 Class * 58 \****************************************************************************************/ 59 60 float EmdL1::getEMDL1(cv::Mat &sig1, cv::Mat &sig2) 61 { 62 // Initialization 63 CV_Assert((sig1.rows==sig2.rows) && (sig1.cols==sig2.cols) && (!sig1.empty()) && (!sig2.empty())); 64 if(!initBaseTrees(sig1.rows, 1)) 65 return -1; 66 67 float *H1=new float[sig1.rows], *H2 = new float[sig2.rows]; 68 for (int ii=0; ii<sig1.rows; ii++) 69 { 70 H1[ii]=sig1.at<float>(ii,0); 71 H2[ii]=sig2.at<float>(ii,0); 72 } 73 74 fillBaseTrees(H1,H2); // Initialize histograms 75 greedySolution(); // Construct an initial Basic Feasible solution 76 initBVTree(); // Initialize BVTree 77 78 // Iteration 79 bool bOptimal = false; 80 m_nItr = 0; 81 while(!bOptimal && m_nItr<nMaxIt) 82 { 83 // Derive U=(u_ij) for row i and column j 84 if(m_nItr==0) updateSubtree(m_pRoot); 85 else updateSubtree(m_pEnter->pChild); 86 87 // Optimality test 88 bOptimal = isOptimal(); 89 90 // Find new solution 91 if(!bOptimal) 92 findNewSolution(); 93 ++m_nItr; 94 } 95 delete [] H1; 96 delete [] H2; 97 // Output the total flow 98 return compuTotalFlow(); 99 } 100 101 void EmdL1::setMaxIteration(int _nMaxIt) 102 { 103 nMaxIt=_nMaxIt; 104 } 105 106 //-- SubFunctions called in the EMD algorithm 107 bool EmdL1::initBaseTrees(int n1, int n2, int n3) 108 { 109 if(binsDim1==n1 && binsDim2==n2 && binsDim3==n3) 110 return true; 111 binsDim1 = n1; 112 binsDim2 = n2; 113 binsDim3 = n3; 114 if(binsDim1==0 || binsDim2==0) dimension = 0; 115 else dimension = (binsDim3==0)?2:3; 116 117 if(dimension==2) 118 { 119 m_Nodes.resize(binsDim1); 120 m_EdgesUp.resize(binsDim1); 121 m_EdgesRight.resize(binsDim1); 122 for(int i1=0; i1<binsDim1; i1++) 123 { 124 m_Nodes[i1].resize(binsDim2); 125 m_EdgesUp[i1].resize(binsDim2); 126 m_EdgesRight[i1].resize(binsDim2); 127 } 128 m_NBVEdges.resize(binsDim1*binsDim2*4+2); 129 m_auxQueue.resize(binsDim1*binsDim2+2); 130 m_fromLoop.resize(binsDim1*binsDim2+2); 131 m_toLoop.resize(binsDim1*binsDim2+2); 132 } 133 else if(dimension==3) 134 { 135 m_3dNodes.resize(binsDim1); 136 m_3dEdgesUp.resize(binsDim1); 137 m_3dEdgesRight.resize(binsDim1); 138 m_3dEdgesDeep.resize(binsDim1); 139 for(int i1=0; i1<binsDim1; i1++) 140 { 141 m_3dNodes[i1].resize(binsDim2); 142 m_3dEdgesUp[i1].resize(binsDim2); 143 m_3dEdgesRight[i1].resize(binsDim2); 144 m_3dEdgesDeep[i1].resize(binsDim2); 145 for(int i2=0; i2<binsDim2; i2++) 146 { 147 m_3dNodes[i1][i2].resize(binsDim3); 148 m_3dEdgesUp[i1][i2].resize(binsDim3); 149 m_3dEdgesRight[i1][i2].resize(binsDim3); 150 m_3dEdgesDeep[i1][i2].resize(binsDim3); 151 } 152 } 153 m_NBVEdges.resize(binsDim1*binsDim2*binsDim3*6+4); 154 m_auxQueue.resize(binsDim1*binsDim2*binsDim3+4); 155 m_fromLoop.resize(binsDim1*binsDim2*binsDim3+4); 156 m_toLoop.resize(binsDim1*binsDim2*binsDim3+2); 157 } 158 else 159 return false; 160 161 return true; 162 } 163 164 bool EmdL1::fillBaseTrees(float *H1, float *H2) 165 { 166 //- Set global counters 167 m_pRoot = NULL; 168 // Graph initialization 169 float *p1 = H1; 170 float *p2 = H2; 171 if(dimension==2) 172 { 173 for(int c=0; c<binsDim2; c++) 174 { 175 for(int r=0; r<binsDim1; r++) 176 { 177 //- initialize nodes and links 178 m_Nodes[r][c].pos[0] = r; 179 m_Nodes[r][c].pos[1] = c; 180 m_Nodes[r][c].d = *(p1++)-*(p2++); 181 m_Nodes[r][c].pParent = NULL; 182 m_Nodes[r][c].pChild = NULL; 183 m_Nodes[r][c].iLevel = -1; 184 185 //- initialize edges 186 // to the right 187 m_EdgesRight[r][c].pParent = &(m_Nodes[r][c]); 188 m_EdgesRight[r][c].pChild = &(m_Nodes[r][(c+1)%binsDim2]); 189 m_EdgesRight[r][c].flow = 0; 190 m_EdgesRight[r][c].iDir = 1; 191 m_EdgesRight[r][c].pNxt = NULL; 192 193 // to the upward 194 m_EdgesUp[r][c].pParent = &(m_Nodes[r][c]); 195 m_EdgesUp[r][c].pChild = &(m_Nodes[(r+1)%binsDim1][c]); 196 m_EdgesUp[r][c].flow = 0; 197 m_EdgesUp[r][c].iDir = 1; 198 m_EdgesUp[r][c].pNxt = NULL; 199 } 200 } 201 } 202 else if(dimension==3) 203 { 204 for(int z=0; z<binsDim3; z++) 205 { 206 for(int c=0; c<binsDim2; c++) 207 { 208 for(int r=0; r<binsDim1; r++) 209 { 210 //- initialize nodes and edges 211 m_3dNodes[r][c][z].pos[0] = r; 212 m_3dNodes[r][c][z].pos[1] = c; 213 m_3dNodes[r][c][z].pos[2] = z; 214 m_3dNodes[r][c][z].d = *(p1++)-*(p2++); 215 m_3dNodes[r][c][z].pParent = NULL; 216 m_3dNodes[r][c][z].pChild = NULL; 217 m_3dNodes[r][c][z].iLevel = -1; 218 219 //- initialize edges 220 // to the upward 221 m_3dEdgesUp[r][c][z].pParent= &(m_3dNodes[r][c][z]); 222 m_3dEdgesUp[r][c][z].pChild = &(m_3dNodes[(r+1)%binsDim1][c][z]); 223 m_3dEdgesUp[r][c][z].flow = 0; 224 m_3dEdgesUp[r][c][z].iDir = 1; 225 m_3dEdgesUp[r][c][z].pNxt = NULL; 226 227 // to the right 228 m_3dEdgesRight[r][c][z].pParent = &(m_3dNodes[r][c][z]); 229 m_3dEdgesRight[r][c][z].pChild = &(m_3dNodes[r][(c+1)%binsDim2][z]); 230 m_3dEdgesRight[r][c][z].flow = 0; 231 m_3dEdgesRight[r][c][z].iDir = 1; 232 m_3dEdgesRight[r][c][z].pNxt = NULL; 233 234 // to the deep 235 m_3dEdgesDeep[r][c][z].pParent = &(m_3dNodes[r][c][z]); 236 m_3dEdgesDeep[r][c][z].pChild = &(m_3dNodes[r][c])[(z+1)%binsDim3]; 237 m_3dEdgesDeep[r][c][z].flow = 0; 238 m_3dEdgesDeep[r][c][z].iDir = 1; 239 m_3dEdgesDeep[r][c][z].pNxt = NULL; 240 } 241 } 242 } 243 } 244 return true; 245 } 246 247 bool EmdL1::greedySolution() 248 { 249 return dimension==2?greedySolution2():greedySolution3(); 250 } 251 252 bool EmdL1::greedySolution2() 253 { 254 //- Prepare auxiliary array, D=H1-H2 255 int c,r; 256 floatArray2D D(binsDim1); 257 for(r=0; r<binsDim1; r++) 258 { 259 D[r].resize(binsDim2); 260 for(c=0; c<binsDim2; c++) D[r][c] = m_Nodes[r][c].d; 261 } 262 // compute integrated values along each dimension 263 std::vector<float> d2s(binsDim2); 264 d2s[0] = 0; 265 for(c=0; c<binsDim2-1; c++) 266 { 267 d2s[c+1] = d2s[c]; 268 for(r=0; r<binsDim1; r++) d2s[c+1]-= D[r][c]; 269 } 270 271 std::vector<float> d1s(binsDim1); 272 d1s[0] = 0; 273 for(r=0; r<binsDim1-1; r++) 274 { 275 d1s[r+1] = d1s[r]; 276 for(c=0; c<binsDim2; c++) d1s[r+1]-= D[r][c]; 277 } 278 279 //- Greedy algorithm for initial solution 280 cvPEmdEdge pBV; 281 float dFlow; 282 bool bUpward = false; 283 nNBV = 0; // number of NON-BV edges 284 285 for(c=0; c<binsDim2-1; c++) 286 for(r=0; r<binsDim1; r++) 287 { 288 dFlow = D[r][c]; 289 bUpward = (r<binsDim1-1) && (fabs(dFlow+d2s[c+1]) > fabs(dFlow+d1s[r+1])); // Move upward or right 290 291 // modify basic variables, record BV and related values 292 if(bUpward) 293 { 294 // move to up 295 pBV = &(m_EdgesUp[r][c]); 296 m_NBVEdges[nNBV++] = &(m_EdgesRight[r][c]); 297 D[r+1][c] += dFlow; // auxilary matrix maintanence 298 d1s[r+1] += dFlow; // auxilary matrix maintanence 299 } 300 else 301 { 302 // move to right, no other choice 303 pBV = &(m_EdgesRight[r][c]); 304 if(r<binsDim1-1) 305 m_NBVEdges[nNBV++] = &(m_EdgesUp[r][c]); 306 307 D[r][c+1] += dFlow; // auxilary matrix maintanence 308 d2s[c+1] += dFlow; // auxilary matrix maintanence 309 } 310 pBV->pParent->pChild = pBV; 311 pBV->flow = fabs(dFlow); 312 pBV->iDir = dFlow>0; // 1:outward, 0:inward 313 } 314 315 //- rightmost column, no choice but move upward 316 c = binsDim2-1; 317 for(r=0; r<binsDim1-1; r++) 318 { 319 dFlow = D[r][c]; 320 pBV = &(m_EdgesUp[r][c]); 321 D[r+1][c] += dFlow; // auxilary matrix maintanence 322 pBV->pParent->pChild= pBV; 323 pBV->flow = fabs(dFlow); 324 pBV->iDir = dFlow>0; // 1:outward, 0:inward 325 } 326 return true; 327 } 328 329 bool EmdL1::greedySolution3() 330 { 331 //- Prepare auxiliary array, D=H1-H2 332 int i1,i2,i3; 333 std::vector<floatArray2D> D(binsDim1); 334 for(i1=0; i1<binsDim1; i1++) 335 { 336 D[i1].resize(binsDim2); 337 for(i2=0; i2<binsDim2; i2++) 338 { 339 D[i1][i2].resize(binsDim3); 340 for(i3=0; i3<binsDim3; i3++) 341 D[i1][i2][i3] = m_3dNodes[i1][i2][i3].d; 342 } 343 } 344 345 // compute integrated values along each dimension 346 std::vector<float> d1s(binsDim1); 347 d1s[0] = 0; 348 for(i1=0; i1<binsDim1-1; i1++) 349 { 350 d1s[i1+1] = d1s[i1]; 351 for(i2=0; i2<binsDim2; i2++) 352 { 353 for(i3=0; i3<binsDim3; i3++) 354 d1s[i1+1] -= D[i1][i2][i3]; 355 } 356 } 357 358 std::vector<float> d2s(binsDim2); 359 d2s[0] = 0; 360 for(i2=0; i2<binsDim2-1; i2++) 361 { 362 d2s[i2+1] = d2s[i2]; 363 for(i1=0; i1<binsDim1; i1++) 364 { 365 for(i3=0; i3<binsDim3; i3++) 366 d2s[i2+1] -= D[i1][i2][i3]; 367 } 368 } 369 370 std::vector<float> d3s(binsDim3); 371 d3s[0] = 0; 372 for(i3=0; i3<binsDim3-1; i3++) 373 { 374 d3s[i3+1] = d3s[i3]; 375 for(i1=0; i1<binsDim1; i1++) 376 { 377 for(i2=0; i2<binsDim2; i2++) 378 d3s[i3+1] -= D[i1][i2][i3]; 379 } 380 } 381 382 //- Greedy algorithm for initial solution 383 cvPEmdEdge pBV; 384 float dFlow, f1,f2,f3; 385 nNBV = 0; // number of NON-BV edges 386 for(i3=0; i3<binsDim3; i3++) 387 { 388 for(i2=0; i2<binsDim2; i2++) 389 { 390 for(i1=0; i1<binsDim1; i1++) 391 { 392 if(i3==binsDim3-1 && i2==binsDim2-1 && i1==binsDim1-1) break; 393 394 //- determine which direction to move, either right or upward 395 dFlow = D[i1][i2][i3]; 396 f1 = (i1<(binsDim1-1))?fabs(dFlow+d1s[i1+1]):std::numeric_limits<float>::max(); 397 f2 = (i2<(binsDim2-1))?fabs(dFlow+d2s[i2+1]):std::numeric_limits<float>::max(); 398 f3 = (i3<(binsDim3-1))?fabs(dFlow+d3s[i3+1]):std::numeric_limits<float>::max(); 399 400 if(f1<f2 && f1<f3) 401 { 402 pBV = &(m_3dEdgesUp[i1][i2][i3]); // up 403 if(i2<binsDim2-1) m_NBVEdges[nNBV++] = &(m_3dEdgesRight[i1][i2][i3]); // right 404 if(i3<binsDim3-1) m_NBVEdges[nNBV++] = &(m_3dEdgesDeep[i1][i2][i3]); // deep 405 D[i1+1][i2][i3] += dFlow; // maintain auxilary matrix 406 d1s[i1+1] += dFlow; 407 } 408 else if(f2<f3) 409 { 410 pBV = &(m_3dEdgesRight[i1][i2][i3]); // right 411 if(i1<binsDim1-1) m_NBVEdges[nNBV++] = &(m_3dEdgesUp[i1][i2][i3]); // up 412 if(i3<binsDim3-1) m_NBVEdges[nNBV++] = &(m_3dEdgesDeep[i1][i2][i3]); // deep 413 D[i1][i2+1][i3] += dFlow; // maintain auxilary matrix 414 d2s[i2+1] += dFlow; 415 } 416 else 417 { 418 pBV = &(m_3dEdgesDeep[i1][i2][i3]); // deep 419 if(i2<binsDim2-1) m_NBVEdges[nNBV++] = &(m_3dEdgesRight[i1][i2][i3]); // right 420 if(i1<binsDim1-1) m_NBVEdges[nNBV++] = &(m_3dEdgesUp[i1][i2][i3]); // up 421 D[i1][i2][i3+1] += dFlow; // maintain auxilary matrix 422 d3s[i3+1] += dFlow; 423 } 424 425 pBV->flow = fabs(dFlow); 426 pBV->iDir = dFlow>0; // 1:outward, 0:inward 427 pBV->pParent->pChild= pBV; 428 } 429 } 430 } 431 return true; 432 } 433 434 void EmdL1::initBVTree() 435 { 436 // initialize BVTree from the initial BF solution 437 //- Using the center of the graph as the root 438 int r = (int)(0.5*binsDim1-.5); 439 int c = (int)(0.5*binsDim2-.5); 440 int z = (int)(0.5*binsDim3-.5); 441 m_pRoot = dimension==2 ? &(m_Nodes[r][c]) : &(m_3dNodes[r][c][z]); 442 m_pRoot->u = 0; 443 m_pRoot->iLevel = 0; 444 m_pRoot->pParent= NULL; 445 m_pRoot->pPEdge = NULL; 446 447 //- Prepare a queue 448 m_auxQueue[0] = m_pRoot; 449 int nQueue = 1; // length of queue 450 int iQHead = 0; // head of queue 451 452 //- Recursively build subtrees 453 cvPEmdEdge pCurE=NULL, pNxtE=NULL; 454 cvPEmdNode pCurN=NULL, pNxtN=NULL; 455 int nBin = binsDim1*binsDim2*std::max(binsDim3,1); 456 while(iQHead<nQueue && nQueue<nBin) 457 { 458 pCurN = m_auxQueue[iQHead++]; // pop out from queue 459 r = pCurN->pos[0]; 460 c = pCurN->pos[1]; 461 z = pCurN->pos[2]; 462 463 // check connection from itself 464 pCurE = pCurN->pChild; // the initial child from initial solution 465 if(pCurE) 466 { 467 pNxtN = pCurE->pChild; 468 pNxtN->pParent = pCurN; 469 pNxtN->pPEdge = pCurE; 470 m_auxQueue[nQueue++] = pNxtN; 471 } 472 473 // check four neighbor nodes 474 int nNB = dimension==2?4:6; 475 for(int k=0;k<nNB;k++) 476 { 477 if(dimension==2) 478 { 479 if(k==0 && c>0) pNxtN = &(m_Nodes[r][c-1]); // left 480 else if(k==1 && r>0) pNxtN = &(m_Nodes[r-1][c]); // down 481 else if(k==2 && c<binsDim2-1) pNxtN = &(m_Nodes[r][c+1]); // right 482 else if(k==3 && r<binsDim1-1) pNxtN = &(m_Nodes[r+1][c]); // up 483 else continue; 484 } 485 else if(dimension==3) 486 { 487 if(k==0 && c>0) pNxtN = &(m_3dNodes[r][c-1][z]); // left 488 else if(k==1 && c<binsDim2-1) pNxtN = &(m_3dNodes[r][c+1][z]); // right 489 else if(k==2 && r>0) pNxtN = &(m_3dNodes[r-1][c][z]); // down 490 else if(k==3 && r<binsDim1-1) pNxtN = &(m_3dNodes[r+1][c][z]); // up 491 else if(k==4 && z>0) pNxtN = &(m_3dNodes[r][c][z-1]); // shallow 492 else if(k==5 && z<binsDim3-1) pNxtN = &(m_3dNodes[r][c][z+1]); // deep 493 else continue; 494 } 495 if(pNxtN != pCurN->pParent) 496 { 497 pNxtE = pNxtN->pChild; 498 if(pNxtE && pNxtE->pChild==pCurN) // has connection 499 { 500 pNxtN->pParent = pCurN; 501 pNxtN->pPEdge = pNxtE; 502 pNxtN->pChild = NULL; 503 m_auxQueue[nQueue++] = pNxtN; 504 505 pNxtE->pParent = pCurN; // reverse direction 506 pNxtE->pChild = pNxtN; 507 pNxtE->iDir = !pNxtE->iDir; 508 509 if(pCurE) pCurE->pNxt = pNxtE; // add to edge list 510 else pCurN->pChild = pNxtE; 511 pCurE = pNxtE; 512 } 513 } 514 } 515 } 516 } 517 518 void EmdL1::updateSubtree(cvPEmdNode pRoot) 519 { 520 // Initialize auxiliary queue 521 m_auxQueue[0] = pRoot; 522 int nQueue = 1; // queue length 523 int iQHead = 0; // head of queue 524 525 // BFS browing 526 cvPEmdNode pCurN=NULL,pNxtN=NULL; 527 cvPEmdEdge pCurE=NULL; 528 while(iQHead<nQueue) 529 { 530 pCurN = m_auxQueue[iQHead++]; // pop out from queue 531 pCurE = pCurN->pChild; 532 533 // browsing all children 534 while(pCurE) 535 { 536 pNxtN = pCurE->pChild; 537 pNxtN->iLevel = pCurN->iLevel+1; 538 pNxtN->u = pCurE->iDir ? (pCurN->u - 1) : (pCurN->u + 1); 539 pCurE = pCurE->pNxt; 540 m_auxQueue[nQueue++] = pNxtN; 541 } 542 } 543 } 544 545 bool EmdL1::isOptimal() 546 { 547 int iC, iMinC = 0; 548 cvPEmdEdge pE; 549 m_pEnter = NULL; 550 m_iEnter = -1; 551 552 // test each NON-BV edges 553 for(int k=0; k<nNBV; ++k) 554 { 555 pE = m_NBVEdges[k]; 556 iC = 1 - pE->pParent->u + pE->pChild->u; 557 if(iC<iMinC) 558 { 559 iMinC = iC; 560 m_iEnter= k; 561 } 562 else 563 { 564 // Try reversing the direction 565 iC = 1 + pE->pParent->u - pE->pChild->u; 566 if(iC<iMinC) 567 { 568 iMinC = iC; 569 m_iEnter= k; 570 } 571 } 572 } 573 574 if(m_iEnter>=0) 575 { 576 m_pEnter = m_NBVEdges[m_iEnter]; 577 if(iMinC == (1 - m_pEnter->pChild->u + m_pEnter->pParent->u)) { 578 // reverse direction 579 cvPEmdNode pN = m_pEnter->pParent; 580 m_pEnter->pParent = m_pEnter->pChild; 581 m_pEnter->pChild = pN; 582 } 583 584 m_pEnter->iDir = 1; 585 } 586 return m_iEnter==-1; 587 } 588 589 void EmdL1::findNewSolution() 590 { 591 // Find loop formed by adding the Enter BV edge. 592 findLoopFromEnterBV(); 593 // Modify flow values along the loop 594 cvPEmdEdge pE = NULL; 595 float minFlow = m_pLeave->flow; 596 int k; 597 for(k=0; k<m_iFrom; k++) 598 { 599 pE = m_fromLoop[k]; 600 if(pE->iDir) pE->flow += minFlow; // outward 601 else pE->flow -= minFlow; // inward 602 } 603 for(k=0; k<m_iTo; k++) 604 { 605 pE = m_toLoop[k]; 606 if(pE->iDir) pE->flow -= minFlow; // outward 607 else pE->flow += minFlow; // inward 608 } 609 610 // Update BV Tree, removing the Leaving-BV edge 611 cvPEmdNode pLParentN = m_pLeave->pParent; 612 cvPEmdNode pLChildN = m_pLeave->pChild; 613 cvPEmdEdge pPreE = pLParentN->pChild; 614 if(pPreE==m_pLeave) 615 { 616 pLParentN->pChild = m_pLeave->pNxt; // Leaving-BV is the first child 617 } 618 else 619 { 620 while(pPreE->pNxt != m_pLeave) 621 pPreE = pPreE->pNxt; 622 pPreE->pNxt = m_pLeave->pNxt; // remove Leaving-BV from child list 623 } 624 pLChildN->pParent = NULL; 625 pLChildN->pPEdge = NULL; 626 627 m_NBVEdges[m_iEnter]= m_pLeave; // put the leaving-BV into the NBV array 628 629 // Add the Enter BV edge 630 cvPEmdNode pEParentN = m_pEnter->pParent; 631 cvPEmdNode pEChildN = m_pEnter->pChild; 632 m_pEnter->flow = minFlow; 633 m_pEnter->pNxt = pEParentN->pChild; // insert the Enter BV as the first child 634 pEParentN->pChild = m_pEnter; // of its parent 635 636 // Recursively update the tree start from pEChildN 637 cvPEmdNode pPreN = pEParentN; 638 cvPEmdNode pCurN = pEChildN; 639 cvPEmdNode pNxtN; 640 cvPEmdEdge pNxtE, pPreE0; 641 pPreE = m_pEnter; 642 while(pCurN) 643 { 644 pNxtN = pCurN->pParent; 645 pNxtE = pCurN->pPEdge; 646 pCurN->pParent = pPreN; 647 pCurN->pPEdge = pPreE; 648 if(pNxtN) 649 { 650 // remove the edge from pNxtN's child list 651 if(pNxtN->pChild==pNxtE) 652 { 653 pNxtN->pChild = pNxtE->pNxt; // first child 654 } 655 else 656 { 657 pPreE0 = pNxtN->pChild; 658 while(pPreE0->pNxt != pNxtE) 659 pPreE0 = pPreE0->pNxt; 660 pPreE0->pNxt = pNxtE->pNxt; // remove Leaving-BV from child list 661 } 662 // reverse the parent-child direction 663 pNxtE->pParent = pCurN; 664 pNxtE->pChild = pNxtN; 665 pNxtE->iDir = !pNxtE->iDir; 666 pNxtE->pNxt = pCurN->pChild; 667 pCurN->pChild = pNxtE; 668 pPreE = pNxtE; 669 pPreN = pCurN; 670 } 671 pCurN = pNxtN; 672 } 673 674 // Update U at the child of the Enter BV 675 pEChildN->u = m_pEnter->iDir?(pEParentN->u-1):(pEParentN->u + 1); 676 pEChildN->iLevel = pEParentN->iLevel+1; 677 } 678 679 void EmdL1::findLoopFromEnterBV() 680 { 681 // Initialize Leaving-BV edge 682 float minFlow = std::numeric_limits<float>::max(); 683 cvPEmdEdge pE = NULL; 684 int iLFlag = 0; // 0: in the FROM list, 1: in the TO list 685 686 // Using two loop list to store the loop nodes 687 cvPEmdNode pFrom = m_pEnter->pParent; 688 cvPEmdNode pTo = m_pEnter->pChild; 689 m_iFrom = 0; 690 m_iTo = 0; 691 m_pLeave = NULL; 692 693 // Trace back to make pFrom and pTo at the same level 694 while(pFrom->iLevel > pTo->iLevel) 695 { 696 pE = pFrom->pPEdge; 697 m_fromLoop[m_iFrom++] = pE; 698 if(!pE->iDir && pE->flow<minFlow) 699 { 700 minFlow = pE->flow; 701 m_pLeave = pE; 702 iLFlag = 0; // 0: in the FROM list 703 } 704 pFrom = pFrom->pParent; 705 } 706 707 while(pTo->iLevel > pFrom->iLevel) 708 { 709 pE = pTo->pPEdge; 710 m_toLoop[m_iTo++] = pE; 711 if(pE->iDir && pE->flow<minFlow) 712 { 713 minFlow = pE->flow; 714 m_pLeave = pE; 715 iLFlag = 1; // 1: in the TO list 716 } 717 pTo = pTo->pParent; 718 } 719 720 // Trace pTo and pFrom simultaneously till find their common ancester 721 while(pTo!=pFrom) 722 { 723 pE = pFrom->pPEdge; 724 m_fromLoop[m_iFrom++] = pE; 725 if(!pE->iDir && pE->flow<minFlow) 726 { 727 minFlow = pE->flow; 728 m_pLeave = pE; 729 iLFlag = 0; // 0: in the FROM list, 1: in the TO list 730 } 731 pFrom = pFrom->pParent; 732 733 pE = pTo->pPEdge; 734 m_toLoop[m_iTo++] = pE; 735 if(pE->iDir && pE->flow<minFlow) 736 { 737 minFlow = pE->flow; 738 m_pLeave = pE; 739 iLFlag = 1; // 0: in the FROM list, 1: in the TO list 740 } 741 pTo = pTo->pParent; 742 } 743 744 // Reverse the direction of the Enter BV edge if necessary 745 if(iLFlag==0) 746 { 747 cvPEmdNode pN = m_pEnter->pParent; 748 m_pEnter->pParent = m_pEnter->pChild; 749 m_pEnter->pChild = pN; 750 m_pEnter->iDir = !m_pEnter->iDir; 751 } 752 } 753 754 float EmdL1::compuTotalFlow() 755 { 756 // Computing the total flow as the final distance 757 float f = 0; 758 759 // Initialize auxiliary queue 760 m_auxQueue[0] = m_pRoot; 761 int nQueue = 1; // length of queue 762 int iQHead = 0; // head of queue 763 764 // BFS browing the tree 765 cvPEmdNode pCurN=NULL,pNxtN=NULL; 766 cvPEmdEdge pCurE=NULL; 767 while(iQHead<nQueue) 768 { 769 pCurN = m_auxQueue[iQHead++]; // pop out from queue 770 pCurE = pCurN->pChild; 771 772 // browsing all children 773 while(pCurE) 774 { 775 f += pCurE->flow; 776 pNxtN = pCurE->pChild; 777 pCurE = pCurE->pNxt; 778 m_auxQueue[nQueue++] = pNxtN; 779 } 780 } 781 return f; 782 } 783 784 /****************************************************************************************\ 785 * EMDL1 Function * 786 \****************************************************************************************/ 787 788 float cv::EMDL1(InputArray _signature1, InputArray _signature2) 789 { 790 Mat signature1 = _signature1.getMat(), signature2 = _signature2.getMat(); 791 EmdL1 emdl1; 792 return emdl1.getEMDL1(signature1, signature2); 793 } 794
