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 
nextEdge(int edge) const46 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 
rotateEdge(int edge,int rotate) const52 int Subdiv2D::rotateEdge(int edge, int rotate) const
53 {
54     return (edge & ~3) + ((edge + rotate) & 3);
55 }
56 
symEdge(int edge) const57 int Subdiv2D::symEdge(int edge) const
58 {
59     return edge ^ 2;
60 }
61 
getEdge(int edge,int nextEdgeType) const62 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 
edgeOrg(int edge,CV_OUT Point2f * orgpt) const69 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 
edgeDst(int edge,CV_OUT Point2f * dstpt) const81 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 
getVertex(int vertex,CV_OUT int * firstEdge) const94 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 
Subdiv2D()103 Subdiv2D::Subdiv2D()
104 {
105     validGeometry = false;
106     freeQEdge = 0;
107     freePoint = 0;
108     recentEdge = 0;
109 }
110 
Subdiv2D(Rect rect)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 
QuadEdge()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 
QuadEdge(int edgeidx)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 
isfree() const139 bool Subdiv2D::QuadEdge::isfree() const
140 {
141     return next[0] <= 0;
142 }
143 
Vertex()144 Subdiv2D::Vertex::Vertex()
145 {
146     firstEdge = 0;
147     type = -1;
148 }
149 
Vertex(Point2f _pt,bool _isvirtual,int _firstEdge)150 Subdiv2D::Vertex::Vertex(Point2f _pt, bool _isvirtual, int _firstEdge)
151 {
152     firstEdge = _firstEdge;
153     type = (int)_isvirtual;
154     pt = _pt;
155 }
156 
isvirtual() const157 bool Subdiv2D::Vertex::isvirtual() const
158 {
159     return type > 0;
160 }
161 
isfree() const162 bool Subdiv2D::Vertex::isfree() const
163 {
164     return type < 0;
165 }
166 
splice(int edgeA,int edgeB)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 
setEdgePoints(int edge,int orgPt,int dstPt)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 
connectEdges(int edgeA,int edgeB)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 
swapEdges(int edge)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 
triangleArea(Point2f a,Point2f b,Point2f c)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 
isRightOf(Point2f pt,int edge) const218 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 
newEdge()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 
deleteEdge(int edge)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 
newPoint(Point2f pt,bool isvirtual,int firstEdge)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 
deletePoint(int vidx)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 
locate(Point2f pt,int & _edge,int & _vertex)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
isPtInCircle3(Point2f pt,Point2f a,Point2f b,Point2f c)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 
insert(Point2f pt)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 
insert(const std::vector<Point2f> & ptvec)480 void Subdiv2D::insert(const std::vector<Point2f>& ptvec)
481 {
482     for( size_t i = 0; i < ptvec.size(); i++ )
483         insert(ptvec[i]);
484 }
485 
initDelaunay(Rect rect)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 
clearVoronoi()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 
computeVoronoiPoint(Point2f org0,Point2f dst0,Point2f org1,Point2f dst1)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 
calcVoronoi()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
isRightOf2(const Point2f & pt,const Point2f & org,const Point2f & diff)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 
findNearest(Point2f pt,Point2f * nearestPt)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 
getEdgeList(std::vector<Vec4f> & edgeList) const709 void Subdiv2D::getEdgeList(std::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 
getTriangleList(std::vector<Vec6f> & triangleList) const726 void Subdiv2D::getTriangleList(std::vector<Vec6f>& triangleList) const
727 {
728     triangleList.clear();
729     int i, total = (int)(qedges.size()*4);
730     std::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 
getVoronoiFacetList(const std::vector<int> & idx,CV_OUT std::vector<std::vector<Point2f>> & facetList,CV_OUT std::vector<Point2f> & facetCenters)750 void Subdiv2D::getVoronoiFacetList(const std::vector<int>& idx,
751                                    CV_OUT std::vector<std::vector<Point2f> >& facetList,
752                                    CV_OUT std::vector<Point2f>& facetCenters)
753 {
754     calcVoronoi();
755     facetList.clear();
756     facetCenters.clear();
757 
758     std::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 
checkSubdiv() const789 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. */
828