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
getEMDL1(cv::Mat & sig1,cv::Mat & sig2)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
setMaxIteration(int _nMaxIt)101 void EmdL1::setMaxIteration(int _nMaxIt)
102 {
103 nMaxIt=_nMaxIt;
104 }
105
106 //-- SubFunctions called in the EMD algorithm
initBaseTrees(int n1,int n2,int n3)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
fillBaseTrees(float * H1,float * H2)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
greedySolution()247 bool EmdL1::greedySolution()
248 {
249 return dimension==2?greedySolution2():greedySolution3();
250 }
251
greedySolution2()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
greedySolution3()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
initBVTree()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
updateSubtree(cvPEmdNode pRoot)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
isOptimal()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
findNewSolution()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
findLoopFromEnterBV()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
compuTotalFlow()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
EMDL1(InputArray _signature1,InputArray _signature2)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