modules/imgproc/src/subdivision2d.cpp

   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 //                        Intel License Agreement
  11 //                For Open Source Computer Vision Library
  12 //
  13 // Copyright (C) 2000, Intel Corporation, all rights reserved.
  14 // Third party copyrights are property of their respective owners.
  15 //
  16 // Redistribution and use in source and binary forms, with or without modification,
  17 // are permitted provided that the following conditions are met:
  18 //
  19 //   * Redistribution's of source code must retain the above copyright notice,
  20 //     this list of conditions and the following disclaimer.
  21 //
  22 //   * Redistribution's in binary form must reproduce the above copyright notice,
  23 //     this list of conditions and the following disclaimer in the documentation
  24 //     and/or other materials provided with the distribution.
  25 //
  26 //   * The name of Intel Corporation may not be used to endorse or promote products
  27 //     derived from this software without specific prior written permission.
  28 //
  29 // This software is provided by the copyright holders and contributors "as is" and
  30 // any express or implied warranties, including, but not limited to, the implied
  31 // warranties of merchantability and fitness for a particular purpose are disclaimed.
  32 // In no event shall the Intel Corporation or contributors be liable for any direct,
  33 // indirect, incidental, special, exemplary, or consequential damages
  34 // (including, but not limited to, procurement of substitute goods or services;
  35 // loss of use, data, or profits; or business interruption) however caused
  36 // and on any theory of liability, whether in contract, strict liability,
  37 // or tort (including negligence or otherwise) arising in any way out of
  38 // the use of this software, even if advised of the possibility of such damage.
  39 //
  40 //M*/
  41 #include "precomp.hpp"
  42
  43 namespace cv
  44 {
  45
  46 int Subdiv2D::nextEdge(int edge) const
  47 {
  48     CV_DbgAssert((size_t)(edge >> 2) < qedges.size());
  49     return qedges[edge >> 2].next[edge & 3];
  50 }
  51
  52 int Subdiv2D::rotateEdge(int edge, int rotate) const
  53 {
  54     return (edge & ~3) + ((edge + rotate) & 3);
  55 }
  56
  57 int Subdiv2D::symEdge(int edge) const
  58 {
  59     return edge ^ 2;
  60 }
  61
  62 int Subdiv2D::getEdge(int edge, int nextEdgeType) const
  63 {
  64     CV_DbgAssert((size_t)(edge >> 2) < qedges.size());
  65     edge = qedges[edge >> 2].next[(edge + nextEdgeType) & 3];
  66     return (edge & ~3) + ((edge + (nextEdgeType >> 4)) & 3);
  67 }
  68
  69 int Subdiv2D::edgeOrg(int edge, CV_OUT Point2f* orgpt) const
  70 {
  71     CV_DbgAssert((size_t)(edge >> 2) < qedges.size());
  72     int vidx = qedges[edge >> 2].pt[edge & 3];
  73     if( orgpt )
  74     {
  75         CV_DbgAssert((size_t)vidx < vtx.size());
  76         *orgpt = vtx[vidx].pt;
  77     }
  78     return vidx;
  79 }
  80
  81 int Subdiv2D::edgeDst(int edge, CV_OUT Point2f* dstpt) const
  82 {
  83     CV_DbgAssert((size_t)(edge >> 2) < qedges.size());
  84     int vidx = qedges[edge >> 2].pt[(edge + 2) & 3];
  85     if( dstpt )
  86     {
  87         CV_DbgAssert((size_t)vidx < vtx.size());
  88         *dstpt = vtx[vidx].pt;
  89     }
  90     return vidx;
  91 }
  92
  93
  94 Point2f Subdiv2D::getVertex(int vertex, CV_OUT int* firstEdge) const
  95 {
  96     CV_DbgAssert((size_t)vertex < vtx.size());
  97     if( firstEdge )
  98         *firstEdge = vtx[vertex].firstEdge;
  99     return vtx[vertex].pt;
 100 }
 101
 102
 103 Subdiv2D::Subdiv2D()
 104 {
 105     validGeometry = false;
 106     freeQEdge = 0;
 107     freePoint = 0;
 108     recentEdge = 0;
 109 }
 110
 111 Subdiv2D::Subdiv2D(Rect rect)
 112 {
 113     validGeometry = false;
 114     freeQEdge = 0;
 115     freePoint = 0;
 116     recentEdge = 0;
 117
 118     initDelaunay(rect);
 119 }
 120
 121
 122 Subdiv2D::QuadEdge::QuadEdge()
 123 {
 124     next[0] = next[1] = next[2] = next[3] = 0;
 125     pt[0] = pt[1] = pt[2] = pt[3] = 0;
 126 }
 127
 128 Subdiv2D::QuadEdge::QuadEdge(int edgeidx)
 129 {
 130     CV_DbgAssert((edgeidx & 3) == 0);
 131     next[0] = edgeidx;
 132     next[1] = edgeidx+3;
 133     next[2] = edgeidx+2;
 134     next[3] = edgeidx+1;
 135
 136     pt[0] = pt[1] = pt[2] = pt[3] = 0;
 137 }
 138
 139 bool Subdiv2D::QuadEdge::isfree() const
 140 {
 141     return next[0] <= 0;
 142 }
 143
 144 Subdiv2D::Vertex::Vertex()
 145 {
 146     firstEdge = 0;
 147     type = -1;
 148 }
 149
 150 Subdiv2D::Vertex::Vertex(Point2f _pt, bool _isvirtual, int _firstEdge)
 151 {
 152     firstEdge = _firstEdge;
 153     type = (int)_isvirtual;
 154     pt = _pt;
 155 }
 156
 157 bool Subdiv2D::Vertex::isvirtual() const
 158 {
 159     return type > 0;
 160 }
 161
 162 bool Subdiv2D::Vertex::isfree() const
 163 {
 164     return type < 0;
 165 }
 166
 167 void Subdiv2D::splice( int edgeA, int edgeB )
 168 {
 169     int& a_next = qedges[edgeA >> 2].next[edgeA & 3];
 170     int& b_next = qedges[edgeB >> 2].next[edgeB & 3];
 171     int a_rot = rotateEdge(a_next, 1);
 172     int b_rot = rotateEdge(b_next, 1);
 173     int& a_rot_next = qedges[a_rot >> 2].next[a_rot & 3];
 174     int& b_rot_next = qedges[b_rot >> 2].next[b_rot & 3];
 175     std::swap(a_next, b_next);
 176     std::swap(a_rot_next, b_rot_next);
 177 }
 178
 179 void Subdiv2D::setEdgePoints(int edge, int orgPt, int dstPt)
 180 {
 181     qedges[edge >> 2].pt[edge & 3] = orgPt;
 182     qedges[edge >> 2].pt[(edge + 2) & 3] = dstPt;
 183     vtx[orgPt].firstEdge = edge;
 184     vtx[dstPt].firstEdge = edge ^ 2;
 185 }
 186
 187 int Subdiv2D::connectEdges( int edgeA, int edgeB )
 188 {
 189     int edge = newEdge();
 190
 191     splice(edge, getEdge(edgeA, NEXT_AROUND_LEFT));
 192     splice(symEdge(edge), edgeB);
 193
 194     setEdgePoints(edge, edgeDst(edgeA), edgeOrg(edgeB));
 195     return edge;
 196 }
 197
 198 void Subdiv2D::swapEdges( int edge )
 199 {
 200     int sedge = symEdge(edge);
 201     int a = getEdge(edge, PREV_AROUND_ORG);
 202     int b = getEdge(sedge, PREV_AROUND_ORG);
 203
 204     splice(edge, a);
 205     splice(sedge, b);
 206
 207     setEdgePoints(edge, edgeDst(a), edgeDst(b));
 208
 209     splice(edge, getEdge(a, NEXT_AROUND_LEFT));
 210     splice(sedge, getEdge(b, NEXT_AROUND_LEFT));
 211 }
 212
 213 static double triangleArea( Point2f a, Point2f b, Point2f c )
 214 {
 215     return ((double)b.x - a.x) * ((double)c.y - a.y) - ((double)b.y - a.y) * ((double)c.x - a.x);
 216 }
 217
 218 int Subdiv2D::isRightOf(Point2f pt, int edge) const
 219 {
 220     Point2f org, dst;
 221     edgeOrg(edge, &org);
 222     edgeDst(edge, &dst);
 223     double cw_area = triangleArea( pt, dst, org );
 224
 225     return (cw_area > 0) - (cw_area < 0);
 226 }
 227
 228 int Subdiv2D::newEdge()
 229 {
 230     if( freeQEdge <= 0 )
 231     {
 232         qedges.push_back(QuadEdge());
 233         freeQEdge = (int)(qedges.size()-1);
 234     }
 235     int edge = freeQEdge*4;
 236     freeQEdge = qedges[edge >> 2].next[1];
 237     qedges[edge >> 2] = QuadEdge(edge);
 238     return edge;
 239 }
 240
 241 void Subdiv2D::deleteEdge(int edge)
 242 {
 243     CV_DbgAssert((size_t)(edge >> 2) < (size_t)qedges.size());
 244     splice( edge, getEdge(edge, PREV_AROUND_ORG) );
 245     int sedge = symEdge(edge);
 246     splice(sedge, getEdge(sedge, PREV_AROUND_ORG) );
 247
 248     edge >>= 2;
 249     qedges[edge].next[0] = 0;
 250     qedges[edge].next[1] = freeQEdge;
 251     freeQEdge = edge;
 252 }
 253
 254 int Subdiv2D::newPoint(Point2f pt, bool isvirtual, int firstEdge)
 255 {
 256     if( freePoint == 0 )
 257     {
 258         vtx.push_back(Vertex());
 259         freePoint = (int)(vtx.size()-1);
 260     }
 261     int vidx = freePoint;
 262     freePoint = vtx[vidx].firstEdge;
 263     vtx[vidx] = Vertex(pt, isvirtual, firstEdge);
 264
 265     return vidx;
 266 }
 267
 268 void Subdiv2D::deletePoint(int vidx)
 269 {
 270     CV_DbgAssert( (size_t)vidx < vtx.size() );
 271     vtx[vidx].firstEdge = freePoint;
 272     vtx[vidx].type = -1;
 273     freePoint = vidx;
 274 }
 275
 276 int Subdiv2D::locate(Point2f pt, int& _edge, int& _vertex)
 277 {
 278     int vertex = 0;
 279
 280     int i, maxEdges = (int)(qedges.size() * 4);
 281
 282     if( qedges.size() < (size_t)4 )
 283         CV_Error( CV_StsError, "Subdivision is empty" );
 284
 285     if( pt.x < topLeft.x || pt.y < topLeft.y || pt.x >= bottomRight.x || pt.y >= bottomRight.y )
 286         CV_Error( CV_StsOutOfRange, "" );
 287
 288     int edge = recentEdge;
 289     CV_Assert(edge > 0);
 290
 291     int location = PTLOC_ERROR;
 292
 293     int right_of_curr = isRightOf(pt, edge);
 294     if( right_of_curr > 0 )
 295     {
 296         edge = symEdge(edge);
 297         right_of_curr = -right_of_curr;
 298     }
 299
 300     for( i = 0; i < maxEdges; i++ )
 301     {
 302         int onext_edge = nextEdge( edge );
 303         int dprev_edge = getEdge( edge, PREV_AROUND_DST );
 304
 305         int right_of_onext = isRightOf( pt, onext_edge );
 306         int right_of_dprev = isRightOf( pt, dprev_edge );
 307
 308         if( right_of_dprev > 0 )
 309         {
 310             if( right_of_onext > 0 || (right_of_onext == 0 && right_of_curr == 0) )
 311             {
 312                 location = PTLOC_INSIDE;
 313                 break;
 314             }
 315             else
 316             {
 317                 right_of_curr = right_of_onext;
 318                 edge = onext_edge;
 319             }
 320         }
 321         else
 322         {
 323             if( right_of_onext > 0 )
 324             {
 325                 if( right_of_dprev == 0 && right_of_curr == 0 )
 326                 {
 327                     location = PTLOC_INSIDE;
 328                     break;
 329                 }
 330                 else
 331                 {
 332                     right_of_curr = right_of_dprev;
 333                     edge = dprev_edge;
 334                 }
 335             }
 336             else if( right_of_curr == 0 &&
 337                     isRightOf( vtx[edgeDst(onext_edge)].pt, edge ) >= 0 )
 338             {
 339                 edge = symEdge( edge );
 340             }
 341             else
 342             {
 343                 right_of_curr = right_of_onext;
 344                 edge = onext_edge;
 345             }
 346         }
 347     }
 348
 349     recentEdge = edge;
 350
 351     if( location == PTLOC_INSIDE )
 352     {
 353         Point2f org_pt, dst_pt;
 354         edgeOrg(edge, &org_pt);
 355         edgeDst(edge, &dst_pt);
 356
 357         double t1 = fabs( pt.x - org_pt.x );
 358         t1 += fabs( pt.y - org_pt.y );
 359         double t2 = fabs( pt.x - dst_pt.x );
 360         t2 += fabs( pt.y - dst_pt.y );
 361         double t3 = fabs( org_pt.x - dst_pt.x );
 362         t3 += fabs( org_pt.y - dst_pt.y );
 363
 364         if( t1 < FLT_EPSILON )
 365         {
 366             location = PTLOC_VERTEX;
 367             vertex = edgeOrg( edge );
 368             edge = 0;
 369         }
 370         else if( t2 < FLT_EPSILON )
 371         {
 372             location = PTLOC_VERTEX;
 373             vertex = edgeDst( edge );
 374             edge = 0;
 375         }
 376         else if( (t1 < t3 || t2 < t3) &&
 377                 fabs( triangleArea( pt, org_pt, dst_pt )) < FLT_EPSILON )
 378         {
 379             location = PTLOC_ON_EDGE;
 380             vertex = 0;
 381         }
 382     }
 383
 384     if( location == PTLOC_ERROR )
 385     {
 386         edge = 0;
 387         vertex = 0;
 388     }
 389
 390     _edge = edge;
 391     _vertex = vertex;
 392
 393     return location;
 394 }
 395
 396
 397 inline int
 398 isPtInCircle3( Point2f pt, Point2f a, Point2f b, Point2f c)
 399 {
 400     const double eps = FLT_EPSILON*0.125;
 401     double val = ((double)a.x * a.x + (double)a.y * a.y) * triangleArea( b, c, pt );
 402     val -= ((double)b.x * b.x + (double)b.y * b.y) * triangleArea( a, c, pt );
 403     val += ((double)c.x * c.x + (double)c.y * c.y) * triangleArea( a, b, pt );
 404     val -= ((double)pt.x * pt.x + (double)pt.y * pt.y) * triangleArea( a, b, c );
 405
 406     return val > eps ? 1 : val < -eps ? -1 : 0;
 407 }
 408
 409
 410 int Subdiv2D::insert(Point2f pt)
 411 {
 412     int curr_point = 0, curr_edge = 0, deleted_edge = 0;
 413     int location = locate( pt, curr_edge, curr_point );
 414
 415     if( location == PTLOC_ERROR )
 416         CV_Error( CV_StsBadSize, "" );
 417
 418     if( location == PTLOC_OUTSIDE_RECT )
 419         CV_Error( CV_StsOutOfRange, "" );
 420
 421     if( location == PTLOC_VERTEX )
 422         return curr_point;
 423
 424     if( location == PTLOC_ON_EDGE )
 425     {
 426         deleted_edge = curr_edge;
 427         recentEdge = curr_edge = getEdge( curr_edge, PREV_AROUND_ORG );
 428         deleteEdge(deleted_edge);
 429     }
 430     else if( location == PTLOC_INSIDE )
 431         ;
 432     else
 433         CV_Error_(CV_StsError, ("Subdiv2D::locate returned invalid location = %d", location) );
 434
 435     assert( curr_edge != 0 );
 436     validGeometry = false;
 437
 438     curr_point = newPoint(pt, false);
 439     int base_edge = newEdge();
 440     int first_point = edgeOrg(curr_edge);
 441     setEdgePoints(base_edge, first_point, curr_point);
 442     splice(base_edge, curr_edge);
 443
 444     do
 445     {
 446         base_edge = connectEdges( curr_edge, symEdge(base_edge) );
 447         curr_edge = getEdge(base_edge, PREV_AROUND_ORG);
 448     }
 449     while( edgeDst(curr_edge) != first_point );
 450
 451     curr_edge = getEdge( base_edge, PREV_AROUND_ORG );
 452
 453     int i, max_edges = (int)(qedges.size()*4);
 454
 455     for( i = 0; i < max_edges; i++ )
 456     {
 457         int temp_dst = 0, curr_org = 0, curr_dst = 0;
 458         int temp_edge = getEdge( curr_edge, PREV_AROUND_ORG );
 459
 460         temp_dst = edgeDst( temp_edge );
 461         curr_org = edgeOrg( curr_edge );
 462         curr_dst = edgeDst( curr_edge );
 463
 464         if( isRightOf( vtx[temp_dst].pt, curr_edge ) > 0 &&
 465            isPtInCircle3( vtx[curr_org].pt, vtx[temp_dst].pt,
 466                          vtx[curr_dst].pt, vtx[curr_point].pt ) < 0 )
 467         {
 468             swapEdges( curr_edge );
 469             curr_edge = getEdge( curr_edge, PREV_AROUND_ORG );
 470         }
 471         else if( curr_org == first_point )
 472             break;
 473         else
 474             curr_edge = getEdge( nextEdge( curr_edge ), PREV_AROUND_LEFT );
 475     }
 476
 477     return curr_point;
 478 }
 479
 480 void Subdiv2D::insert(const vector<Point2f>& ptvec)
 481 {
 482     for( size_t i = 0; i < ptvec.size(); i++ )
 483         insert(ptvec[i]);
 484 }
 485
 486 void Subdiv2D::initDelaunay( Rect rect )
 487 {
 488     float big_coord = 3.f * MAX( rect.width, rect.height );
 489     float rx = (float)rect.x;
 490     float ry = (float)rect.y;
 491
 492     vtx.clear();
 493     qedges.clear();
 494
 495     recentEdge = 0;
 496     validGeometry = false;
 497
 498     topLeft = Point2f( rx, ry );
 499     bottomRight = Point2f( rx + rect.width, ry + rect.height );
 500
 501     Point2f ppA( rx + big_coord, ry );
 502     Point2f ppB( rx, ry + big_coord );
 503     Point2f ppC( rx - big_coord, ry - big_coord );
 504
 505     vtx.push_back(Vertex());
 506     qedges.push_back(QuadEdge());
 507
 508     freeQEdge = 0;
 509     freePoint = 0;
 510
 511     int pA = newPoint(ppA, false);
 512     int pB = newPoint(ppB, false);
 513     int pC = newPoint(ppC, false);
 514
 515     int edge_AB = newEdge();
 516     int edge_BC = newEdge();
 517     int edge_CA = newEdge();
 518
 519     setEdgePoints( edge_AB, pA, pB );
 520     setEdgePoints( edge_BC, pB, pC );
 521     setEdgePoints( edge_CA, pC, pA );
 522
 523     splice( edge_AB, symEdge( edge_CA ));
 524     splice( edge_BC, symEdge( edge_AB ));
 525     splice( edge_CA, symEdge( edge_BC ));
 526
 527     recentEdge = edge_AB;
 528 }
 529
 530
 531 void Subdiv2D::clearVoronoi()
 532 {
 533     size_t i, total = qedges.size();
 534
 535     for( i = 0; i < total; i++ )
 536         qedges[i].pt[1] = qedges[i].pt[3] = 0;
 537
 538     total = vtx.size();
 539     for( i = 0; i < total; i++ )
 540     {
 541         if( vtx[i].isvirtual() )
 542             deletePoint((int)i);
 543     }
 544
 545     validGeometry = false;
 546 }
 547
 548
 549 static Point2f computeVoronoiPoint(Point2f org0, Point2f dst0, Point2f org1, Point2f dst1)
 550 {
 551     double a0 = dst0.x - org0.x;
 552     double b0 = dst0.y - org0.y;
 553     double c0 = -0.5*(a0 * (dst0.x + org0.x) + b0 * (dst0.y + org0.y));
 554
 555     double a1 = dst1.x - org1.x;
 556     double b1 = dst1.y - org1.y;
 557     double c1 = -0.5*(a1 * (dst1.x + org1.x) + b1 * (dst1.y + org1.y));
 558
 559     double det = a0 * b1 - a1 * b0;
 560
 561     if( det != 0 )
 562     {
 563         det = 1. / det;
 564         return Point2f((float) ((b0 * c1 - b1 * c0) * det),
 565                        (float) ((a1 * c0 - a0 * c1) * det));
 566     }
 567
 568     return Point2f(FLT_MAX, FLT_MAX);
 569 }
 570
 571
 572 void Subdiv2D::calcVoronoi()
 573 {
 574     // check if it is already calculated
 575     if( validGeometry )
 576         return;
 577
 578     clearVoronoi();
 579     int i, total = (int)qedges.size();
 580
 581     // loop through all quad-edges, except for the first 3 (#1, #2, #3 - 0 is reserved for "NULL" pointer)
 582     for( i = 4; i < total; i++ )
 583     {
 584         QuadEdge& quadedge = qedges[i];
 585
 586         if( quadedge.isfree() )
 587             continue;
 588
 589         int edge0 = (int)(i*4);
 590         Point2f org0, dst0, org1, dst1;
 591
 592         if( !quadedge.pt[3] )
 593         {
 594             int edge1 = getEdge( edge0, NEXT_AROUND_LEFT );
 595             int edge2 = getEdge( edge1, NEXT_AROUND_LEFT );
 596
 597             edgeOrg(edge0, &org0);
 598             edgeDst(edge0, &dst0);
 599             edgeOrg(edge1, &org1);
 600             edgeDst(edge1, &dst1);
 601
 602             Point2f virt_point = computeVoronoiPoint(org0, dst0, org1, dst1);
 603
 604             if( fabs( virt_point.x ) < FLT_MAX * 0.5 &&
 605                fabs( virt_point.y ) < FLT_MAX * 0.5 )
 606             {
 607                 quadedge.pt[3] = qedges[edge1 >> 2].pt[3 - (edge1 & 2)] =
 608                 qedges[edge2 >> 2].pt[3 - (edge2 & 2)] = newPoint(virt_point, true);
 609             }
 610         }
 611
 612         if( !quadedge.pt[1] )
 613         {
 614             int edge1 = getEdge( edge0, NEXT_AROUND_RIGHT );
 615             int edge2 = getEdge( edge1, NEXT_AROUND_RIGHT );
 616
 617             edgeOrg(edge0, &org0);
 618             edgeDst(edge0, &dst0);
 619             edgeOrg(edge1, &org1);
 620             edgeDst(edge1, &dst1);
 621
 622             Point2f virt_point = computeVoronoiPoint(org0, dst0, org1, dst1);
 623
 624             if( fabs( virt_point.x ) < FLT_MAX * 0.5 &&
 625                fabs( virt_point.y ) < FLT_MAX * 0.5 )
 626             {
 627                 quadedge.pt[1] = qedges[edge1 >> 2].pt[1 + (edge1 & 2)] =
 628                 qedges[edge2 >> 2].pt[1 + (edge2 & 2)] = newPoint(virt_point, true);
 629             }
 630         }
 631     }
 632
 633     validGeometry = true;
 634 }
 635
 636
 637 static int
 638 isRightOf2( const Point2f& pt, const Point2f& org, const Point2f& diff )
 639 {
 640     double cw_area = ((double)org.x - pt.x)*diff.y - ((double)org.y - pt.y)*diff.x;
 641     return (cw_area > 0) - (cw_area < 0);
 642 }
 643
 644
 645 int Subdiv2D::findNearest(Point2f pt, Point2f* nearestPt)
 646 {
 647     if( !validGeometry )
 648         calcVoronoi();
 649
 650     int vertex = 0, edge = 0;
 651     int loc = locate( pt, edge, vertex );
 652
 653     if( loc != PTLOC_ON_EDGE && loc != PTLOC_INSIDE )
 654         return vertex;
 655
 656     vertex = 0;
 657
 658     Point2f start;
 659     edgeOrg(edge, &start);
 660     Point2f diff = pt - start;
 661
 662     edge = rotateEdge(edge, 1);
 663
 664     int i, total = (int)vtx.size();
 665
 666     for( i = 0; i < total; i++ )
 667     {
 668         Point2f t;
 669
 670         for(;;)
 671         {
 672             CV_Assert( edgeDst(edge, &t) > 0 );
 673             if( isRightOf2( t, start, diff ) >= 0 )
 674                 break;
 675
 676             edge = getEdge( edge, NEXT_AROUND_LEFT );
 677         }
 678
 679         for(;;)
 680         {
 681             CV_Assert( edgeOrg( edge, &t ) > 0 );
 682
 683             if( isRightOf2( t, start, diff ) < 0 )
 684                 break;
 685
 686             edge = getEdge( edge, PREV_AROUND_LEFT );
 687         }
 688
 689         Point2f tempDiff;
 690         edgeDst(edge, &tempDiff);
 691         edgeOrg(edge, &t);
 692         tempDiff -= t;
 693
 694         if( isRightOf2( pt, t, tempDiff ) >= 0 )
 695         {
 696             vertex = edgeOrg(rotateEdge( edge, 3 ));
 697             break;
 698         }
 699
 700         edge = symEdge( edge );
 701     }
 702
 703     if( nearestPt && vertex > 0 )
 704         *nearestPt = vtx[vertex].pt;
 705
 706     return vertex;
 707 }
 708
 709 void Subdiv2D::getEdgeList(vector<Vec4f>& edgeList) const
 710 {
 711     edgeList.clear();
 712
 713     for( size_t i = 4; i < qedges.size(); i++ )
 714     {
 715         if( qedges[i].isfree() )
 716             continue;
 717         if( qedges[i].pt[0] > 0 && qedges[i].pt[2] > 0 )
 718         {
 719             Point2f org = vtx[qedges[i].pt[0]].pt;
 720             Point2f dst = vtx[qedges[i].pt[2]].pt;
 721             edgeList.push_back(Vec4f(org.x, org.y, dst.x, dst.y));
 722         }
 723     }
 724 }
 725
 726 void Subdiv2D::getTriangleList(vector<Vec6f>& triangleList) const
 727 {
 728     triangleList.clear();
 729     int i, total = (int)(qedges.size()*4);
 730     vector<bool> edgemask(total, false);
 731
 732     for( i = 4; i < total; i += 2 )
 733     {
 734         if( edgemask[i] )
 735             continue;
 736         Point2f a, b, c;
 737         int edge = i;
 738         edgeOrg(edge, &a);
 739         edgemask[edge] = true;
 740         edge = getEdge(edge, NEXT_AROUND_LEFT);
 741         edgeOrg(edge, &b);
 742         edgemask[edge] = true;
 743         edge = getEdge(edge, NEXT_AROUND_LEFT);
 744         edgeOrg(edge, &c);
 745         edgemask[edge] = true;
 746         triangleList.push_back(Vec6f(a.x, a.y, b.x, b.y, c.x, c.y));
 747     }
 748 }
 749
 750 void Subdiv2D::getVoronoiFacetList(const vector<int>& idx,
 751                                    CV_OUT vector<vector<Point2f> >& facetList,
 752                                    CV_OUT vector<Point2f>& facetCenters)
 753 {
 754     calcVoronoi();
 755     facetList.clear();
 756     facetCenters.clear();
 757
 758     vector<Point2f> buf;
 759
 760     size_t i, total;
 761     if( idx.empty() )
 762         i = 4, total = vtx.size();
 763     else
 764         i = 0, total = idx.size();
 765
 766     for( ; i < total; i++ )
 767     {
 768         int k = idx.empty() ? (int)i : idx[i];
 769
 770         if( vtx[k].isfree() || vtx[k].isvirtual() )
 771             continue;
 772         int edge = rotateEdge(vtx[k].firstEdge, 1), t = edge;
 773
 774         // gather points
 775         buf.clear();
 776         do
 777         {
 778             buf.push_back(vtx[edgeOrg(t)].pt);
 779             t = getEdge( t, NEXT_AROUND_LEFT );
 780         }
 781         while( t != edge );
 782
 783         facetList.push_back(buf);
 784         facetCenters.push_back(vtx[k].pt);
 785     }
 786 }
 787
 788
 789 void Subdiv2D::checkSubdiv() const
 790 {
 791     int i, j, total = (int)qedges.size();
 792
 793     for( i = 0; i < total; i++ )
 794     {
 795         const QuadEdge& qe = qedges[i];
 796
 797         if( qe.isfree() )
 798             continue;
 799
 800         for( j = 0; j < 4; j++ )
 801         {
 802             int e = (int)(i*4 + j);
 803             int o_next = nextEdge(e);
 804             int o_prev = getEdge(e, PREV_AROUND_ORG );
 805             int d_prev = getEdge(e, PREV_AROUND_DST );
 806             int d_next = getEdge(e, NEXT_AROUND_DST );
 807
 808             // check points
 809             CV_Assert( edgeOrg(e) == edgeOrg(o_next));
 810             CV_Assert( edgeOrg(e) == edgeOrg(o_prev));
 811             CV_Assert( edgeDst(e) == edgeDst(d_next));
 812             CV_Assert( edgeDst(e) == edgeDst(d_prev));
 813
 814             if( j % 2 == 0 )
 815             {
 816                 CV_Assert( edgeDst(o_next) == edgeOrg(d_prev));
 817                 CV_Assert( edgeDst(o_prev) == edgeOrg(d_next));
 818                 CV_Assert( getEdge(getEdge(getEdge(e,NEXT_AROUND_LEFT),NEXT_AROUND_LEFT),NEXT_AROUND_LEFT) == e );
 819                 CV_Assert( getEdge(getEdge(getEdge(e,NEXT_AROUND_RIGHT),NEXT_AROUND_RIGHT),NEXT_AROUND_RIGHT) == e);
 820             }
 821         }
 822     }
 823 }
 824
 825 }
 826
 827 /* End of file. */