1 //Copyright (C) 2011 by Ivan Fratric
2 //
3 //Permission is hereby granted, free of charge, to any person obtaining a copy
4 //of this software and associated documentation files (the "Software"), to deal
5 //in the Software without restriction, including without limitation the rights
6 //to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
7 //copies of the Software, and to permit persons to whom the Software is
8 //furnished to do so, subject to the following conditions:
9 //
10 //The above copyright notice and this permission notice shall be included in
11 //all copies or substantial portions of the Software.
12 //
13 //THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
14 //IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
15 //FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
16 //AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
17 //LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
18 //OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
19 //THE SOFTWARE.
20
21
22 #include <stdio.h>
23 #include <string.h>
24 #include <math.h>
25 #include <list>
26 #include <algorithm>
27 #include <set>
28 #include <vector>
29 #include <stdexcept>
30
31 using namespace std;
32
33 #include "polypartition.h"
34
35 #define TPPL_VERTEXTYPE_REGULAR 0
36 #define TPPL_VERTEXTYPE_START 1
37 #define TPPL_VERTEXTYPE_END 2
38 #define TPPL_VERTEXTYPE_SPLIT 3
39 #define TPPL_VERTEXTYPE_MERGE 4
40
TPPLPoly()41 TPPLPoly::TPPLPoly() {
42 hole = false;
43 numpoints = 0;
44 points = NULL;
45 }
46
~TPPLPoly()47 TPPLPoly::~TPPLPoly() {
48 if(points) delete [] points;
49 }
50
Clear()51 void TPPLPoly::Clear() {
52 if(points) delete [] points;
53 hole = false;
54 numpoints = 0;
55 points = NULL;
56 }
57
Init(long numpoints)58 void TPPLPoly::Init(long numpoints) {
59 Clear();
60 this->numpoints = numpoints;
61 points = new TPPLPoint[numpoints];
62 }
63
Triangle(TPPLPoint & p1,TPPLPoint & p2,TPPLPoint & p3)64 void TPPLPoly::Triangle(TPPLPoint &p1, TPPLPoint &p2, TPPLPoint &p3) {
65 Init(3);
66 points[0] = p1;
67 points[1] = p2;
68 points[2] = p3;
69 }
70
TPPLPoly(const TPPLPoly & src)71 TPPLPoly::TPPLPoly(const TPPLPoly &src) : TPPLPoly() {
72 hole = src.hole;
73 numpoints = src.numpoints;
74
75 if(numpoints > 0) {
76 points = new TPPLPoint[numpoints];
77 memcpy(points, src.points, numpoints*sizeof(TPPLPoint));
78 }
79 }
80
operator =(const TPPLPoly & src)81 TPPLPoly& TPPLPoly::operator=(const TPPLPoly &src) {
82 Clear();
83 hole = src.hole;
84 numpoints = src.numpoints;
85
86 if(numpoints > 0) {
87 points = new TPPLPoint[numpoints];
88 memcpy(points, src.points, numpoints*sizeof(TPPLPoint));
89 }
90
91 return *this;
92 }
93
GetOrientation() const94 int TPPLPoly::GetOrientation() const {
95 long i1,i2;
96 tppl_float area = 0;
97 for(i1=0; i1<numpoints; i1++) {
98 i2 = i1+1;
99 if(i2 == numpoints) i2 = 0;
100 area += points[i1].x * points[i2].y - points[i1].y * points[i2].x;
101 }
102 if(area>0) return TPPL_CCW;
103 if(area<0) return TPPL_CW;
104 return 0;
105 }
106
SetOrientation(int orientation)107 void TPPLPoly::SetOrientation(int orientation) {
108 int polyorientation = GetOrientation();
109 if(polyorientation&&(polyorientation!=orientation)) {
110 Invert();
111 }
112 }
113
Invert()114 void TPPLPoly::Invert() {
115 std::reverse(points, points + numpoints);
116 }
117
PartitionVertex()118 TPPLPartition::PartitionVertex::PartitionVertex() : previous(NULL), next(NULL) {
119
120 }
121
Normalize(const TPPLPoint & p)122 TPPLPoint TPPLPartition::Normalize(const TPPLPoint &p) {
123 TPPLPoint r;
124 tppl_float n = sqrt(p.x*p.x + p.y*p.y);
125 if(n!=0) {
126 r = p/n;
127 } else {
128 r.x = 0;
129 r.y = 0;
130 }
131 return r;
132 }
133
Distance(const TPPLPoint & p1,const TPPLPoint & p2)134 tppl_float TPPLPartition::Distance(const TPPLPoint &p1, const TPPLPoint &p2) {
135 tppl_float dx,dy;
136 dx = p2.x - p1.x;
137 dy = p2.y - p1.y;
138 return(sqrt(dx*dx + dy*dy));
139 }
140
141 //checks if two lines intersect
Intersects(TPPLPoint & p11,TPPLPoint & p12,TPPLPoint & p21,TPPLPoint & p22)142 int TPPLPartition::Intersects(TPPLPoint &p11, TPPLPoint &p12, TPPLPoint &p21, TPPLPoint &p22) {
143 if((p11.x == p21.x)&&(p11.y == p21.y)) return 0;
144 if((p11.x == p22.x)&&(p11.y == p22.y)) return 0;
145 if((p12.x == p21.x)&&(p12.y == p21.y)) return 0;
146 if((p12.x == p22.x)&&(p12.y == p22.y)) return 0;
147
148 TPPLPoint v1ort,v2ort,v;
149 tppl_float dot11,dot12,dot21,dot22;
150
151 v1ort.x = p12.y-p11.y;
152 v1ort.y = p11.x-p12.x;
153
154 v2ort.x = p22.y-p21.y;
155 v2ort.y = p21.x-p22.x;
156
157 v = p21-p11;
158 dot21 = v.x*v1ort.x + v.y*v1ort.y;
159 v = p22-p11;
160 dot22 = v.x*v1ort.x + v.y*v1ort.y;
161
162 v = p11-p21;
163 dot11 = v.x*v2ort.x + v.y*v2ort.y;
164 v = p12-p21;
165 dot12 = v.x*v2ort.x + v.y*v2ort.y;
166
167 if(dot11*dot12>0) return 0;
168 if(dot21*dot22>0) return 0;
169
170 return 1;
171 }
172
173 //removes holes from inpolys by merging them with non-holes
RemoveHoles(TPPLPolyList * inpolys,TPPLPolyList * outpolys)174 int TPPLPartition::RemoveHoles(TPPLPolyList *inpolys, TPPLPolyList *outpolys) {
175 TPPLPolyList polys;
176 TPPLPolyList::iterator holeiter,polyiter,iter,iter2;
177 long i,i2,holepointindex,polypointindex;
178 TPPLPoint holepoint,polypoint,bestpolypoint;
179 TPPLPoint linep1,linep2;
180 TPPLPoint v1,v2;
181 TPPLPoly newpoly;
182 bool hasholes;
183 bool pointvisible;
184 bool pointfound;
185
186 //check for trivial case (no holes)
187 hasholes = false;
188 for(iter = inpolys->begin(); iter!=inpolys->end(); iter++) {
189 if(iter->IsHole()) {
190 hasholes = true;
191 break;
192 }
193 }
194 if(!hasholes) {
195 for(iter = inpolys->begin(); iter!=inpolys->end(); iter++) {
196 outpolys->push_back(*iter);
197 }
198 return 1;
199 }
200
201 polys = *inpolys;
202
203 while(1) {
204 //find the hole point with the largest x
205 hasholes = false;
206 for(iter = polys.begin(); iter!=polys.end(); iter++) {
207 if(!iter->IsHole()) continue;
208
209 if(!hasholes) {
210 hasholes = true;
211 holeiter = iter;
212 holepointindex = 0;
213 }
214
215 for(i=0; i < iter->GetNumPoints(); i++) {
216 if(iter->GetPoint(i).x > holeiter->GetPoint(holepointindex).x) {
217 holeiter = iter;
218 holepointindex = i;
219 }
220 }
221 }
222 if(!hasholes) break;
223 holepoint = holeiter->GetPoint(holepointindex);
224
225 pointfound = false;
226 for(iter = polys.begin(); iter!=polys.end(); iter++) {
227 if(iter->IsHole()) continue;
228 for(i=0; i < iter->GetNumPoints(); i++) {
229 if(iter->GetPoint(i).x <= holepoint.x) continue;
230 if(!InCone(iter->GetPoint((i+iter->GetNumPoints()-1)%(iter->GetNumPoints())),
231 iter->GetPoint(i),
232 iter->GetPoint((i+1)%(iter->GetNumPoints())),
233 holepoint))
234 continue;
235 polypoint = iter->GetPoint(i);
236 if(pointfound) {
237 v1 = Normalize(polypoint-holepoint);
238 v2 = Normalize(bestpolypoint-holepoint);
239 if(v2.x > v1.x) continue;
240 }
241 pointvisible = true;
242 for(iter2 = polys.begin(); iter2!=polys.end(); iter2++) {
243 if(iter2->IsHole()) continue;
244 for(i2=0; i2 < iter2->GetNumPoints(); i2++) {
245 linep1 = iter2->GetPoint(i2);
246 linep2 = iter2->GetPoint((i2+1)%(iter2->GetNumPoints()));
247 if(Intersects(holepoint,polypoint,linep1,linep2)) {
248 pointvisible = false;
249 break;
250 }
251 }
252 if(!pointvisible) break;
253 }
254 if(pointvisible) {
255 pointfound = true;
256 bestpolypoint = polypoint;
257 polyiter = iter;
258 polypointindex = i;
259 }
260 }
261 }
262
263 if(!pointfound) return 0;
264
265 newpoly.Init(holeiter->GetNumPoints() + polyiter->GetNumPoints() + 2);
266 i2 = 0;
267 for(i=0;i<=polypointindex;i++) {
268 newpoly[i2] = polyiter->GetPoint(i);
269 i2++;
270 }
271 for(i=0;i<=holeiter->GetNumPoints();i++) {
272 newpoly[i2] = holeiter->GetPoint((i+holepointindex)%holeiter->GetNumPoints());
273 i2++;
274 }
275 for(i=polypointindex;i<polyiter->GetNumPoints();i++) {
276 newpoly[i2] = polyiter->GetPoint(i);
277 i2++;
278 }
279
280 polys.erase(holeiter);
281 polys.erase(polyiter);
282 polys.push_back(newpoly);
283 }
284
285 for(iter = polys.begin(); iter!=polys.end(); iter++) {
286 outpolys->push_back(*iter);
287 }
288
289 return 1;
290 }
291
IsConvex(TPPLPoint & p1,TPPLPoint & p2,TPPLPoint & p3)292 bool TPPLPartition::IsConvex(TPPLPoint& p1, TPPLPoint& p2, TPPLPoint& p3) {
293 tppl_float tmp;
294 tmp = (p3.y-p1.y)*(p2.x-p1.x)-(p3.x-p1.x)*(p2.y-p1.y);
295 if(tmp>0) return 1;
296 else return 0;
297 }
298
IsReflex(TPPLPoint & p1,TPPLPoint & p2,TPPLPoint & p3)299 bool TPPLPartition::IsReflex(TPPLPoint& p1, TPPLPoint& p2, TPPLPoint& p3) {
300 tppl_float tmp;
301 tmp = (p3.y-p1.y)*(p2.x-p1.x)-(p3.x-p1.x)*(p2.y-p1.y);
302 if(tmp<0) return 1;
303 else return 0;
304 }
305
IsInside(TPPLPoint & p1,TPPLPoint & p2,TPPLPoint & p3,TPPLPoint & p)306 bool TPPLPartition::IsInside(TPPLPoint& p1, TPPLPoint& p2, TPPLPoint& p3, TPPLPoint &p) {
307 if(IsConvex(p1,p,p2)) return false;
308 if(IsConvex(p2,p,p3)) return false;
309 if(IsConvex(p3,p,p1)) return false;
310 return true;
311 }
312
InCone(TPPLPoint & p1,TPPLPoint & p2,TPPLPoint & p3,TPPLPoint & p)313 bool TPPLPartition::InCone(TPPLPoint &p1, TPPLPoint &p2, TPPLPoint &p3, TPPLPoint &p) {
314 bool convex;
315
316 convex = IsConvex(p1,p2,p3);
317
318 if(convex) {
319 if(!IsConvex(p1,p2,p)) return false;
320 if(!IsConvex(p2,p3,p)) return false;
321 return true;
322 } else {
323 if(IsConvex(p1,p2,p)) return true;
324 if(IsConvex(p2,p3,p)) return true;
325 return false;
326 }
327 }
328
InCone(PartitionVertex * v,TPPLPoint & p)329 bool TPPLPartition::InCone(PartitionVertex *v, TPPLPoint &p) {
330 TPPLPoint p1,p2,p3;
331
332 p1 = v->previous->p;
333 p2 = v->p;
334 p3 = v->next->p;
335
336 return InCone(p1,p2,p3,p);
337 }
338
UpdateVertexReflexity(PartitionVertex * v)339 void TPPLPartition::UpdateVertexReflexity(PartitionVertex *v) {
340 PartitionVertex *v1 = NULL,*v3 = NULL;
341 v1 = v->previous;
342 v3 = v->next;
343 v->isConvex = !IsReflex(v1->p,v->p,v3->p);
344 }
345
UpdateVertex(PartitionVertex * v,PartitionVertex * vertices,long numvertices)346 void TPPLPartition::UpdateVertex(PartitionVertex *v, PartitionVertex *vertices, long numvertices) {
347 long i;
348 PartitionVertex *v1 = NULL,*v3 = NULL;
349 TPPLPoint vec1,vec3;
350
351 v1 = v->previous;
352 v3 = v->next;
353
354 v->isConvex = IsConvex(v1->p,v->p,v3->p);
355
356 vec1 = Normalize(v1->p - v->p);
357 vec3 = Normalize(v3->p - v->p);
358 v->angle = vec1.x*vec3.x + vec1.y*vec3.y;
359
360 if(v->isConvex) {
361 v->isEar = true;
362 for(i=0;i<numvertices;i++) {
363 if((vertices[i].p.x==v->p.x)&&(vertices[i].p.y==v->p.y)) continue;
364 if((vertices[i].p.x==v1->p.x)&&(vertices[i].p.y==v1->p.y)) continue;
365 if((vertices[i].p.x==v3->p.x)&&(vertices[i].p.y==v3->p.y)) continue;
366 if(IsInside(v1->p,v->p,v3->p,vertices[i].p)) {
367 v->isEar = false;
368 break;
369 }
370 }
371 } else {
372 v->isEar = false;
373 }
374 }
375
376 //triangulation by ear removal
Triangulate_EC(TPPLPoly * poly,TPPLPolyList * triangles)377 int TPPLPartition::Triangulate_EC(TPPLPoly *poly, TPPLPolyList *triangles) {
378 if(!poly->Valid()) return 0;
379
380 long numvertices;
381 PartitionVertex *vertices = NULL;
382 PartitionVertex *ear = NULL;
383 TPPLPoly triangle;
384 long i,j;
385 bool earfound;
386
387 if(poly->GetNumPoints() < 3) return 0;
388 if(poly->GetNumPoints() == 3) {
389 triangles->push_back(*poly);
390 return 1;
391 }
392
393 numvertices = poly->GetNumPoints();
394
395 vertices = new PartitionVertex[numvertices];
396 for(i=0;i<numvertices;i++) {
397 vertices[i].isActive = true;
398 vertices[i].p = poly->GetPoint(i);
399 if(i==(numvertices-1)) vertices[i].next=&(vertices[0]);
400 else vertices[i].next=&(vertices[i+1]);
401 if(i==0) vertices[i].previous = &(vertices[numvertices-1]);
402 else vertices[i].previous = &(vertices[i-1]);
403 }
404 for(i=0;i<numvertices;i++) {
405 UpdateVertex(&vertices[i],vertices,numvertices);
406 }
407
408 for(i=0;i<numvertices-3;i++) {
409 earfound = false;
410 //find the most extruded ear
411 for(j=0;j<numvertices;j++) {
412 if(!vertices[j].isActive) continue;
413 if(!vertices[j].isEar) continue;
414 if(!earfound) {
415 earfound = true;
416 ear = &(vertices[j]);
417 } else {
418 if(vertices[j].angle > ear->angle) {
419 ear = &(vertices[j]);
420 }
421 }
422 }
423 if(!earfound) {
424 delete [] vertices;
425 return 0;
426 }
427
428 triangle.Triangle(ear->previous->p,ear->p,ear->next->p);
429 triangles->push_back(triangle);
430
431 ear->isActive = false;
432 ear->previous->next = ear->next;
433 ear->next->previous = ear->previous;
434
435 if(i==numvertices-4) break;
436
437 UpdateVertex(ear->previous,vertices,numvertices);
438 UpdateVertex(ear->next,vertices,numvertices);
439 }
440 for(i=0;i<numvertices;i++) {
441 if(vertices[i].isActive) {
442 triangle.Triangle(vertices[i].previous->p,vertices[i].p,vertices[i].next->p);
443 triangles->push_back(triangle);
444 break;
445 }
446 }
447
448 delete [] vertices;
449
450 return 1;
451 }
452
Triangulate_EC(TPPLPolyList * inpolys,TPPLPolyList * triangles)453 int TPPLPartition::Triangulate_EC(TPPLPolyList *inpolys, TPPLPolyList *triangles) {
454 TPPLPolyList outpolys;
455 TPPLPolyList::iterator iter;
456
457 if(!RemoveHoles(inpolys,&outpolys)) return 0;
458 for(iter=outpolys.begin();iter!=outpolys.end();iter++) {
459 if(!Triangulate_EC(&(*iter),triangles)) return 0;
460 }
461 return 1;
462 }
463
ConvexPartition_HM(TPPLPoly * poly,TPPLPolyList * parts)464 int TPPLPartition::ConvexPartition_HM(TPPLPoly *poly, TPPLPolyList *parts) {
465 if(!poly->Valid()) return 0;
466
467 TPPLPolyList triangles;
468 TPPLPolyList::iterator iter1,iter2;
469 TPPLPoly *poly1 = NULL,*poly2 = NULL;
470 TPPLPoly newpoly;
471 TPPLPoint d1,d2,p1,p2,p3;
472 long i11,i12,i21,i22,i13,i23,j,k;
473 bool isdiagonal;
474 long numreflex;
475
476 //check if the poly is already convex
477 numreflex = 0;
478 for(i11=0;i11<poly->GetNumPoints();i11++) {
479 if(i11==0) i12 = poly->GetNumPoints()-1;
480 else i12=i11-1;
481 if(i11==(poly->GetNumPoints()-1)) i13=0;
482 else i13=i11+1;
483 if(IsReflex(poly->GetPoint(i12),poly->GetPoint(i11),poly->GetPoint(i13))) {
484 numreflex = 1;
485 break;
486 }
487 }
488 if(numreflex == 0) {
489 parts->push_back(*poly);
490 return 1;
491 }
492
493 if(!Triangulate_EC(poly,&triangles)) return 0;
494
495 for(iter1 = triangles.begin(); iter1 != triangles.end(); iter1++) {
496 poly1 = &(*iter1);
497 for(i11=0;i11<poly1->GetNumPoints();i11++) {
498 d1 = poly1->GetPoint(i11);
499 i12 = (i11+1)%(poly1->GetNumPoints());
500 d2 = poly1->GetPoint(i12);
501
502 isdiagonal = false;
503 for(iter2 = iter1; iter2 != triangles.end(); iter2++) {
504 if(iter1 == iter2) continue;
505 poly2 = &(*iter2);
506
507 for(i21=0;i21<poly2->GetNumPoints();i21++) {
508 if((d2.x != poly2->GetPoint(i21).x)||(d2.y != poly2->GetPoint(i21).y)) continue;
509 i22 = (i21+1)%(poly2->GetNumPoints());
510 if((d1.x != poly2->GetPoint(i22).x)||(d1.y != poly2->GetPoint(i22).y)) continue;
511 isdiagonal = true;
512 break;
513 }
514 if(isdiagonal) break;
515 }
516
517 if(!isdiagonal) continue;
518
519 p2 = poly1->GetPoint(i11);
520 if(i11 == 0) i13 = poly1->GetNumPoints()-1;
521 else i13 = i11-1;
522 p1 = poly1->GetPoint(i13);
523 if(i22 == (poly2->GetNumPoints()-1)) i23 = 0;
524 else i23 = i22+1;
525 p3 = poly2->GetPoint(i23);
526
527 if(!IsConvex(p1,p2,p3)) continue;
528
529 p2 = poly1->GetPoint(i12);
530 if(i12 == (poly1->GetNumPoints()-1)) i13 = 0;
531 else i13 = i12+1;
532 p3 = poly1->GetPoint(i13);
533 if(i21 == 0) i23 = poly2->GetNumPoints()-1;
534 else i23 = i21-1;
535 p1 = poly2->GetPoint(i23);
536
537 if(!IsConvex(p1,p2,p3)) continue;
538
539 newpoly.Init(poly1->GetNumPoints()+poly2->GetNumPoints()-2);
540 k = 0;
541 for(j=i12;j!=i11;j=(j+1)%(poly1->GetNumPoints())) {
542 newpoly[k] = poly1->GetPoint(j);
543 k++;
544 }
545 for(j=i22;j!=i21;j=(j+1)%(poly2->GetNumPoints())) {
546 newpoly[k] = poly2->GetPoint(j);
547 k++;
548 }
549
550 triangles.erase(iter2);
551 *iter1 = newpoly;
552 poly1 = &(*iter1);
553 i11 = -1;
554
555 continue;
556 }
557 }
558
559 for(iter1 = triangles.begin(); iter1 != triangles.end(); iter1++) {
560 parts->push_back(*iter1);
561 }
562
563 return 1;
564 }
565
ConvexPartition_HM(TPPLPolyList * inpolys,TPPLPolyList * parts)566 int TPPLPartition::ConvexPartition_HM(TPPLPolyList *inpolys, TPPLPolyList *parts) {
567 TPPLPolyList outpolys;
568 TPPLPolyList::iterator iter;
569
570 if(!RemoveHoles(inpolys,&outpolys)) return 0;
571 for(iter=outpolys.begin();iter!=outpolys.end();iter++) {
572 if(!ConvexPartition_HM(&(*iter),parts)) return 0;
573 }
574 return 1;
575 }
576
577 //minimum-weight polygon triangulation by dynamic programming
578 //O(n^3) time complexity
579 //O(n^2) space complexity
Triangulate_OPT(TPPLPoly * poly,TPPLPolyList * triangles)580 int TPPLPartition::Triangulate_OPT(TPPLPoly *poly, TPPLPolyList *triangles) {
581 if(!poly->Valid()) return 0;
582
583 long i,j,k,gap,n;
584 DPState **dpstates = NULL;
585 TPPLPoint p1,p2,p3,p4;
586 long bestvertex;
587 tppl_float weight,minweight,d1,d2;
588 Diagonal diagonal,newdiagonal;
589 DiagonalList diagonals;
590 TPPLPoly triangle;
591 int ret = 1;
592
593 n = poly->GetNumPoints();
594 dpstates = new DPState *[n];
595 for(i=1;i<n;i++) {
596 dpstates[i] = new DPState[i];
597 }
598
599 //init states and visibility
600 for(i=0;i<(n-1);i++) {
601 p1 = poly->GetPoint(i);
602 for(j=i+1;j<n;j++) {
603 dpstates[j][i].visible = true;
604 dpstates[j][i].weight = 0;
605 dpstates[j][i].bestvertex = -1;
606 if(j!=(i+1)) {
607 p2 = poly->GetPoint(j);
608
609 //visibility check
610 if(i==0) p3 = poly->GetPoint(n-1);
611 else p3 = poly->GetPoint(i-1);
612 if(i==(n-1)) p4 = poly->GetPoint(0);
613 else p4 = poly->GetPoint(i+1);
614 if(!InCone(p3,p1,p4,p2)) {
615 dpstates[j][i].visible = false;
616 continue;
617 }
618
619 if(j==0) p3 = poly->GetPoint(n-1);
620 else p3 = poly->GetPoint(j-1);
621 if(j==(n-1)) p4 = poly->GetPoint(0);
622 else p4 = poly->GetPoint(j+1);
623 if(!InCone(p3,p2,p4,p1)) {
624 dpstates[j][i].visible = false;
625 continue;
626 }
627
628 for(k=0;k<n;k++) {
629 p3 = poly->GetPoint(k);
630 if(k==(n-1)) p4 = poly->GetPoint(0);
631 else p4 = poly->GetPoint(k+1);
632 if(Intersects(p1,p2,p3,p4)) {
633 dpstates[j][i].visible = false;
634 break;
635 }
636 }
637 }
638 }
639 }
640 dpstates[n-1][0].visible = true;
641 dpstates[n-1][0].weight = 0;
642 dpstates[n-1][0].bestvertex = -1;
643
644 for(gap = 2; gap<n; gap++) {
645 for(i=0; i<(n-gap); i++) {
646 j = i+gap;
647 if(!dpstates[j][i].visible) continue;
648 bestvertex = -1;
649 for(k=(i+1);k<j;k++) {
650 if(!dpstates[k][i].visible) continue;
651 if(!dpstates[j][k].visible) continue;
652
653 if(k<=(i+1)) d1=0;
654 else d1 = Distance(poly->GetPoint(i),poly->GetPoint(k));
655 if(j<=(k+1)) d2=0;
656 else d2 = Distance(poly->GetPoint(k),poly->GetPoint(j));
657
658 weight = dpstates[k][i].weight + dpstates[j][k].weight + d1 + d2;
659
660 if((bestvertex == -1)||(weight<minweight)) {
661 bestvertex = k;
662 minweight = weight;
663 }
664 }
665 if(bestvertex == -1) {
666 for(i=1;i<n;i++) {
667 delete [] dpstates[i];
668 }
669 delete [] dpstates;
670
671 return 0;
672 }
673
674 dpstates[j][i].bestvertex = bestvertex;
675 dpstates[j][i].weight = minweight;
676 }
677 }
678
679 newdiagonal.index1 = 0;
680 newdiagonal.index2 = n-1;
681 diagonals.push_back(newdiagonal);
682 while(!diagonals.empty()) {
683 diagonal = *(diagonals.begin());
684 diagonals.pop_front();
685 bestvertex = dpstates[diagonal.index2][diagonal.index1].bestvertex;
686 if(bestvertex == -1) {
687 ret = 0;
688 break;
689 }
690 triangle.Triangle(poly->GetPoint(diagonal.index1),poly->GetPoint(bestvertex),poly->GetPoint(diagonal.index2));
691 triangles->push_back(triangle);
692 if(bestvertex > (diagonal.index1+1)) {
693 newdiagonal.index1 = diagonal.index1;
694 newdiagonal.index2 = bestvertex;
695 diagonals.push_back(newdiagonal);
696 }
697 if(diagonal.index2 > (bestvertex+1)) {
698 newdiagonal.index1 = bestvertex;
699 newdiagonal.index2 = diagonal.index2;
700 diagonals.push_back(newdiagonal);
701 }
702 }
703
704 for(i=1;i<n;i++) {
705 delete [] dpstates[i];
706 }
707 delete [] dpstates;
708
709 return ret;
710 }
711
UpdateState(long a,long b,long w,long i,long j,DPState2 ** dpstates)712 void TPPLPartition::UpdateState(long a, long b, long w, long i, long j, DPState2 **dpstates) {
713 Diagonal newdiagonal;
714 DiagonalList *pairs = NULL;
715 long w2;
716
717 w2 = dpstates[a][b].weight;
718 if(w>w2) return;
719
720 pairs = &(dpstates[a][b].pairs);
721 newdiagonal.index1 = i;
722 newdiagonal.index2 = j;
723
724 if(w<w2) {
725 pairs->clear();
726 pairs->push_front(newdiagonal);
727 dpstates[a][b].weight = w;
728 } else {
729 if((!pairs->empty())&&(i <= pairs->begin()->index1)) return;
730 while((!pairs->empty())&&(pairs->begin()->index2 >= j)) pairs->pop_front();
731 pairs->push_front(newdiagonal);
732 }
733 }
734
TypeA(long i,long j,long k,PartitionVertex * vertices,DPState2 ** dpstates)735 void TPPLPartition::TypeA(long i, long j, long k, PartitionVertex *vertices, DPState2 **dpstates) {
736 DiagonalList *pairs = NULL;
737 DiagonalList::iterator iter,lastiter;
738 long top;
739 long w;
740
741 if(!dpstates[i][j].visible) return;
742 top = j;
743 w = dpstates[i][j].weight;
744 if(k-j > 1) {
745 if (!dpstates[j][k].visible) return;
746 w += dpstates[j][k].weight + 1;
747 }
748 if(j-i > 1) {
749 pairs = &(dpstates[i][j].pairs);
750 iter = pairs->end();
751 lastiter = pairs->end();
752 while(iter!=pairs->begin()) {
753 iter--;
754 if(!IsReflex(vertices[iter->index2].p,vertices[j].p,vertices[k].p)) lastiter = iter;
755 else break;
756 }
757 if(lastiter == pairs->end()) w++;
758 else {
759 if(IsReflex(vertices[k].p,vertices[i].p,vertices[lastiter->index1].p)) w++;
760 else top = lastiter->index1;
761 }
762 }
763 UpdateState(i,k,w,top,j,dpstates);
764 }
765
TypeB(long i,long j,long k,PartitionVertex * vertices,DPState2 ** dpstates)766 void TPPLPartition::TypeB(long i, long j, long k, PartitionVertex *vertices, DPState2 **dpstates) {
767 DiagonalList *pairs = NULL;
768 DiagonalList::iterator iter,lastiter;
769 long top;
770 long w;
771
772 if(!dpstates[j][k].visible) return;
773 top = j;
774 w = dpstates[j][k].weight;
775
776 if (j-i > 1) {
777 if (!dpstates[i][j].visible) return;
778 w += dpstates[i][j].weight + 1;
779 }
780 if (k-j > 1) {
781 pairs = &(dpstates[j][k].pairs);
782
783 iter = pairs->begin();
784 if((!pairs->empty())&&(!IsReflex(vertices[i].p,vertices[j].p,vertices[iter->index1].p))) {
785 lastiter = iter;
786 while(iter!=pairs->end()) {
787 if(!IsReflex(vertices[i].p,vertices[j].p,vertices[iter->index1].p)) {
788 lastiter = iter;
789 iter++;
790 }
791 else break;
792 }
793 if(IsReflex(vertices[lastiter->index2].p,vertices[k].p,vertices[i].p)) w++;
794 else top = lastiter->index2;
795 } else w++;
796 }
797 UpdateState(i,k,w,j,top,dpstates);
798 }
799
ConvexPartition_OPT(TPPLPoly * poly,TPPLPolyList * parts)800 int TPPLPartition::ConvexPartition_OPT(TPPLPoly *poly, TPPLPolyList *parts) {
801 if(!poly->Valid()) return 0;
802
803 TPPLPoint p1,p2,p3,p4;
804 PartitionVertex *vertices = NULL;
805 DPState2 **dpstates = NULL;
806 long i,j,k,n,gap;
807 DiagonalList diagonals,diagonals2;
808 Diagonal diagonal,newdiagonal;
809 DiagonalList *pairs = NULL,*pairs2 = NULL;
810 DiagonalList::iterator iter,iter2;
811 int ret;
812 TPPLPoly newpoly;
813 vector<long> indices;
814 vector<long>::iterator iiter;
815 bool ijreal,jkreal;
816
817 n = poly->GetNumPoints();
818 vertices = new PartitionVertex[n];
819
820 dpstates = new DPState2 *[n];
821 for(i=0;i<n;i++) {
822 dpstates[i] = new DPState2[n];
823 }
824
825 //init vertex information
826 for(i=0;i<n;i++) {
827 vertices[i].p = poly->GetPoint(i);
828 vertices[i].isActive = true;
829 if(i==0) vertices[i].previous = &(vertices[n-1]);
830 else vertices[i].previous = &(vertices[i-1]);
831 if(i==(poly->GetNumPoints()-1)) vertices[i].next = &(vertices[0]);
832 else vertices[i].next = &(vertices[i+1]);
833 }
834 for(i=1;i<n;i++) {
835 UpdateVertexReflexity(&(vertices[i]));
836 }
837
838 //init states and visibility
839 for(i=0;i<(n-1);i++) {
840 p1 = poly->GetPoint(i);
841 for(j=i+1;j<n;j++) {
842 dpstates[i][j].visible = true;
843 if(j==i+1) {
844 dpstates[i][j].weight = 0;
845 } else {
846 dpstates[i][j].weight = 2147483647;
847 }
848 if(j!=(i+1)) {
849 p2 = poly->GetPoint(j);
850
851 //visibility check
852 if(!InCone(&vertices[i],p2)) {
853 dpstates[i][j].visible = false;
854 continue;
855 }
856 if(!InCone(&vertices[j],p1)) {
857 dpstates[i][j].visible = false;
858 continue;
859 }
860
861 for(k=0;k<n;k++) {
862 p3 = poly->GetPoint(k);
863 if(k==(n-1)) p4 = poly->GetPoint(0);
864 else p4 = poly->GetPoint(k+1);
865 if(Intersects(p1,p2,p3,p4)) {
866 dpstates[i][j].visible = false;
867 break;
868 }
869 }
870 }
871 }
872 }
873 for(i=0;i<(n-2);i++) {
874 j = i+2;
875 if(dpstates[i][j].visible) {
876 dpstates[i][j].weight = 0;
877 newdiagonal.index1 = i+1;
878 newdiagonal.index2 = i+1;
879 dpstates[i][j].pairs.push_back(newdiagonal);
880 }
881 }
882
883 dpstates[0][n-1].visible = true;
884 vertices[0].isConvex = false; //by convention
885
886 for(gap=3; gap<n; gap++) {
887 for(i=0;i<n-gap;i++) {
888 if(vertices[i].isConvex) continue;
889 k = i+gap;
890 if(dpstates[i][k].visible) {
891 if(!vertices[k].isConvex) {
892 for(j=i+1;j<k;j++) TypeA(i,j,k,vertices,dpstates);
893 } else {
894 for(j=i+1;j<(k-1);j++) {
895 if(vertices[j].isConvex) continue;
896 TypeA(i,j,k,vertices,dpstates);
897 }
898 TypeA(i,k-1,k,vertices,dpstates);
899 }
900 }
901 }
902 for(k=gap;k<n;k++) {
903 if(vertices[k].isConvex) continue;
904 i = k-gap;
905 if((vertices[i].isConvex)&&(dpstates[i][k].visible)) {
906 TypeB(i,i+1,k,vertices,dpstates);
907 for(j=i+2;j<k;j++) {
908 if(vertices[j].isConvex) continue;
909 TypeB(i,j,k,vertices,dpstates);
910 }
911 }
912 }
913 }
914
915
916 //recover solution
917 ret = 1;
918 newdiagonal.index1 = 0;
919 newdiagonal.index2 = n-1;
920 diagonals.push_front(newdiagonal);
921 while(!diagonals.empty()) {
922 diagonal = *(diagonals.begin());
923 diagonals.pop_front();
924 if((diagonal.index2 - diagonal.index1) <=1) continue;
925 pairs = &(dpstates[diagonal.index1][diagonal.index2].pairs);
926 if(pairs->empty()) {
927 ret = 0;
928 break;
929 }
930 if(!vertices[diagonal.index1].isConvex) {
931 iter = pairs->end();
932 iter--;
933 j = iter->index2;
934 newdiagonal.index1 = j;
935 newdiagonal.index2 = diagonal.index2;
936 diagonals.push_front(newdiagonal);
937 if((j - diagonal.index1)>1) {
938 if(iter->index1 != iter->index2) {
939 pairs2 = &(dpstates[diagonal.index1][j].pairs);
940 while(1) {
941 if(pairs2->empty()) {
942 ret = 0;
943 break;
944 }
945 iter2 = pairs2->end();
946 iter2--;
947 if(iter->index1 != iter2->index1) pairs2->pop_back();
948 else break;
949 }
950 if(ret == 0) break;
951 }
952 newdiagonal.index1 = diagonal.index1;
953 newdiagonal.index2 = j;
954 diagonals.push_front(newdiagonal);
955 }
956 } else {
957 iter = pairs->begin();
958 j = iter->index1;
959 newdiagonal.index1 = diagonal.index1;
960 newdiagonal.index2 = j;
961 diagonals.push_front(newdiagonal);
962 if((diagonal.index2 - j) > 1) {
963 if(iter->index1 != iter->index2) {
964 pairs2 = &(dpstates[j][diagonal.index2].pairs);
965 while(1) {
966 if(pairs2->empty()) {
967 ret = 0;
968 break;
969 }
970 iter2 = pairs2->begin();
971 if(iter->index2 != iter2->index2) pairs2->pop_front();
972 else break;
973 }
974 if(ret == 0) break;
975 }
976 newdiagonal.index1 = j;
977 newdiagonal.index2 = diagonal.index2;
978 diagonals.push_front(newdiagonal);
979 }
980 }
981 }
982
983 if(ret == 0) {
984 for(i=0;i<n;i++) {
985 delete [] dpstates[i];
986 }
987 delete [] dpstates;
988 delete [] vertices;
989
990 return ret;
991 }
992
993 newdiagonal.index1 = 0;
994 newdiagonal.index2 = n-1;
995 diagonals.push_front(newdiagonal);
996 while(!diagonals.empty()) {
997 diagonal = *(diagonals.begin());
998 diagonals.pop_front();
999 if((diagonal.index2 - diagonal.index1) <= 1) continue;
1000
1001 indices.clear();
1002 diagonals2.clear();
1003 indices.push_back(diagonal.index1);
1004 indices.push_back(diagonal.index2);
1005 diagonals2.push_front(diagonal);
1006
1007 while(!diagonals2.empty()) {
1008 diagonal = *(diagonals2.begin());
1009 diagonals2.pop_front();
1010 if((diagonal.index2 - diagonal.index1) <= 1) continue;
1011 ijreal = true;
1012 jkreal = true;
1013 pairs = &(dpstates[diagonal.index1][diagonal.index2].pairs);
1014 if(!vertices[diagonal.index1].isConvex) {
1015 iter = pairs->end();
1016 iter--;
1017 j = iter->index2;
1018 if(iter->index1 != iter->index2) ijreal = false;
1019 } else {
1020 iter = pairs->begin();
1021 j = iter->index1;
1022 if(iter->index1 != iter->index2) jkreal = false;
1023 }
1024
1025 newdiagonal.index1 = diagonal.index1;
1026 newdiagonal.index2 = j;
1027 if(ijreal) {
1028 diagonals.push_back(newdiagonal);
1029 } else {
1030 diagonals2.push_back(newdiagonal);
1031 }
1032
1033 newdiagonal.index1 = j;
1034 newdiagonal.index2 = diagonal.index2;
1035 if(jkreal) {
1036 diagonals.push_back(newdiagonal);
1037 } else {
1038 diagonals2.push_back(newdiagonal);
1039 }
1040
1041 indices.push_back(j);
1042 }
1043
1044 std::sort(indices.begin(), indices.end());
1045 newpoly.Init((long)indices.size());
1046 k=0;
1047 for(iiter = indices.begin();iiter!=indices.end();iiter++) {
1048 newpoly[k] = vertices[*iiter].p;
1049 k++;
1050 }
1051 parts->push_back(newpoly);
1052 }
1053
1054 for(i=0;i<n;i++) {
1055 delete [] dpstates[i];
1056 }
1057 delete [] dpstates;
1058 delete [] vertices;
1059
1060 return ret;
1061 }
1062
1063 //triangulates a set of polygons by first partitioning them into monotone polygons
1064 //O(n*log(n)) time complexity, O(n) space complexity
1065 //the algorithm used here is outlined in the book
1066 //"Computational Geometry: Algorithms and Applications"
1067 //by Mark de Berg, Otfried Cheong, Marc van Kreveld and Mark Overmars
MonotonePartition(TPPLPolyList * inpolys,TPPLPolyList * monotonePolys)1068 int TPPLPartition::MonotonePartition(TPPLPolyList *inpolys, TPPLPolyList *monotonePolys) {
1069 TPPLPolyList::iterator iter;
1070 MonotoneVertex *vertices = NULL;
1071 long i,numvertices,vindex,vindex2,newnumvertices,maxnumvertices;
1072 long polystartindex, polyendindex;
1073 TPPLPoly *poly = NULL;
1074 MonotoneVertex *v = NULL,*v2 = NULL,*vprev = NULL,*vnext = NULL;
1075 ScanLineEdge newedge;
1076 bool error = false;
1077
1078 numvertices = 0;
1079 for(iter = inpolys->begin(); iter != inpolys->end(); iter++) {
1080 if(!iter->Valid()) return 0;
1081 numvertices += iter->GetNumPoints();
1082 }
1083
1084 maxnumvertices = numvertices*3;
1085 vertices = new MonotoneVertex[maxnumvertices];
1086 newnumvertices = numvertices;
1087
1088 polystartindex = 0;
1089 for(iter = inpolys->begin(); iter != inpolys->end(); iter++) {
1090 poly = &(*iter);
1091 polyendindex = polystartindex + poly->GetNumPoints()-1;
1092 for(i=0;i<poly->GetNumPoints();i++) {
1093 vertices[i+polystartindex].p = poly->GetPoint(i);
1094 if(i==0) vertices[i+polystartindex].previous = polyendindex;
1095 else vertices[i+polystartindex].previous = i+polystartindex-1;
1096 if(i==(poly->GetNumPoints()-1)) vertices[i+polystartindex].next = polystartindex;
1097 else vertices[i+polystartindex].next = i+polystartindex+1;
1098 }
1099 polystartindex = polyendindex+1;
1100 }
1101
1102 //construct the priority queue
1103 long *priority = new long [numvertices];
1104 for(i=0;i<numvertices;i++) priority[i] = i;
1105 std::sort(priority,&(priority[numvertices]),VertexSorter(vertices));
1106
1107 //determine vertex types
1108 char *vertextypes = new char[maxnumvertices];
1109 for(i=0;i<numvertices;i++) {
1110 v = &(vertices[i]);
1111 vprev = &(vertices[v->previous]);
1112 vnext = &(vertices[v->next]);
1113
1114 if(Below(vprev->p,v->p)&&Below(vnext->p,v->p)) {
1115 if(IsConvex(vnext->p,vprev->p,v->p)) {
1116 vertextypes[i] = TPPL_VERTEXTYPE_START;
1117 } else {
1118 vertextypes[i] = TPPL_VERTEXTYPE_SPLIT;
1119 }
1120 } else if(Below(v->p,vprev->p)&&Below(v->p,vnext->p)) {
1121 if(IsConvex(vnext->p,vprev->p,v->p))
1122 {
1123 vertextypes[i] = TPPL_VERTEXTYPE_END;
1124 } else {
1125 vertextypes[i] = TPPL_VERTEXTYPE_MERGE;
1126 }
1127 } else {
1128 vertextypes[i] = TPPL_VERTEXTYPE_REGULAR;
1129 }
1130 }
1131
1132 //helpers
1133 long *helpers = new long[maxnumvertices];
1134
1135 //binary search tree that holds edges intersecting the scanline
1136 //note that while set doesn't actually have to be implemented as a tree
1137 //complexity requirements for operations are the same as for the balanced binary search tree
1138 set<ScanLineEdge> edgeTree;
1139 //store iterators to the edge tree elements
1140 //this makes deleting existing edges much faster
1141 set<ScanLineEdge>::iterator *edgeTreeIterators,edgeIter;
1142 edgeTreeIterators = new set<ScanLineEdge>::iterator[maxnumvertices];
1143 pair<set<ScanLineEdge>::iterator,bool> edgeTreeRet;
1144 for(i = 0; i<numvertices; i++) edgeTreeIterators[i] = edgeTree.end();
1145
1146 //for each vertex
1147 for(i=0;i<numvertices;i++) {
1148 vindex = priority[i];
1149 v = &(vertices[vindex]);
1150 vindex2 = vindex;
1151 v2 = v;
1152
1153 //depending on the vertex type, do the appropriate action
1154 //comments in the following sections are copied from "Computational Geometry: Algorithms and Applications"
1155 switch(vertextypes[vindex]) {
1156 case TPPL_VERTEXTYPE_START:
1157 //Insert ei in T and set helper(ei) to vi.
1158 newedge.p1 = v->p;
1159 newedge.p2 = vertices[v->next].p;
1160 newedge.index = vindex;
1161 edgeTreeRet = edgeTree.insert(newedge);
1162 edgeTreeIterators[vindex] = edgeTreeRet.first;
1163 helpers[vindex] = vindex;
1164 break;
1165
1166 case TPPL_VERTEXTYPE_END:
1167 if (edgeTreeIterators[v->previous] == edgeTree.end()) {
1168 error = true;
1169 break;
1170 }
1171 //if helper(ei-1) is a merge vertex
1172 if(vertextypes[helpers[v->previous]]==TPPL_VERTEXTYPE_MERGE) {
1173 //Insert the diagonal connecting vi to helper(ei-1) in D.
1174 AddDiagonal(vertices,&newnumvertices,vindex,helpers[v->previous],
1175 vertextypes, edgeTreeIterators, &edgeTree, helpers);
1176 }
1177 //Delete ei-1 from T
1178 edgeTree.erase(edgeTreeIterators[v->previous]);
1179 break;
1180
1181 case TPPL_VERTEXTYPE_SPLIT:
1182 //Search in T to find the edge e j directly left of vi.
1183 newedge.p1 = v->p;
1184 newedge.p2 = v->p;
1185 edgeIter = edgeTree.lower_bound(newedge);
1186 if(edgeIter == edgeTree.begin()) {
1187 error = true;
1188 break;
1189 }
1190 edgeIter--;
1191 //Insert the diagonal connecting vi to helper(ej) in D.
1192 AddDiagonal(vertices,&newnumvertices,vindex,helpers[edgeIter->index],
1193 vertextypes, edgeTreeIterators, &edgeTree, helpers);
1194 vindex2 = newnumvertices-2;
1195 v2 = &(vertices[vindex2]);
1196 //helper(e j)�vi
1197 helpers[edgeIter->index] = vindex;
1198 //Insert ei in T and set helper(ei) to vi.
1199 newedge.p1 = v2->p;
1200 newedge.p2 = vertices[v2->next].p;
1201 newedge.index = vindex2;
1202 edgeTreeRet = edgeTree.insert(newedge);
1203 edgeTreeIterators[vindex2] = edgeTreeRet.first;
1204 helpers[vindex2] = vindex2;
1205 break;
1206
1207 case TPPL_VERTEXTYPE_MERGE:
1208 if (edgeTreeIterators[v->previous] == edgeTree.end()) {
1209 error = true;
1210 break;
1211 }
1212 //if helper(ei-1) is a merge vertex
1213 if(vertextypes[helpers[v->previous]]==TPPL_VERTEXTYPE_MERGE) {
1214 //Insert the diagonal connecting vi to helper(ei-1) in D.
1215 AddDiagonal(vertices,&newnumvertices,vindex,helpers[v->previous],
1216 vertextypes, edgeTreeIterators, &edgeTree, helpers);
1217 vindex2 = newnumvertices-2;
1218 v2 = &(vertices[vindex2]);
1219 }
1220 //Delete ei-1 from T.
1221 edgeTree.erase(edgeTreeIterators[v->previous]);
1222 //Search in T to find the edge e j directly left of vi.
1223 newedge.p1 = v->p;
1224 newedge.p2 = v->p;
1225 edgeIter = edgeTree.lower_bound(newedge);
1226 if(edgeIter == edgeTree.begin()) {
1227 error = true;
1228 break;
1229 }
1230 edgeIter--;
1231 //if helper(ej) is a merge vertex
1232 if(vertextypes[helpers[edgeIter->index]]==TPPL_VERTEXTYPE_MERGE) {
1233 //Insert the diagonal connecting vi to helper(e j) in D.
1234 AddDiagonal(vertices,&newnumvertices,vindex2,helpers[edgeIter->index],
1235 vertextypes, edgeTreeIterators, &edgeTree, helpers);
1236 }
1237 //helper(e j)�vi
1238 helpers[edgeIter->index] = vindex2;
1239 break;
1240
1241 case TPPL_VERTEXTYPE_REGULAR:
1242 //if the interior of P lies to the right of vi
1243 if(Below(v->p,vertices[v->previous].p)) {
1244 if (edgeTreeIterators[v->previous] == edgeTree.end()) {
1245 error = true;
1246 break;
1247 }
1248 //if helper(ei-1) is a merge vertex
1249 if(vertextypes[helpers[v->previous]]==TPPL_VERTEXTYPE_MERGE) {
1250 //Insert the diagonal connecting vi to helper(ei-1) in D.
1251 AddDiagonal(vertices,&newnumvertices,vindex,helpers[v->previous],
1252 vertextypes, edgeTreeIterators, &edgeTree, helpers);
1253 vindex2 = newnumvertices-2;
1254 v2 = &(vertices[vindex2]);
1255 }
1256 //Delete ei-1 from T.
1257 edgeTree.erase(edgeTreeIterators[v->previous]);
1258 //Insert ei in T and set helper(ei) to vi.
1259 newedge.p1 = v2->p;
1260 newedge.p2 = vertices[v2->next].p;
1261 newedge.index = vindex2;
1262 edgeTreeRet = edgeTree.insert(newedge);
1263 edgeTreeIterators[vindex2] = edgeTreeRet.first;
1264 helpers[vindex2] = vindex;
1265 } else {
1266 //Search in T to find the edge ej directly left of vi.
1267 newedge.p1 = v->p;
1268 newedge.p2 = v->p;
1269 edgeIter = edgeTree.lower_bound(newedge);
1270 if(edgeIter == edgeTree.begin()) {
1271 error = true;
1272 break;
1273 }
1274 edgeIter--;
1275 //if helper(ej) is a merge vertex
1276 if(vertextypes[helpers[edgeIter->index]]==TPPL_VERTEXTYPE_MERGE) {
1277 //Insert the diagonal connecting vi to helper(e j) in D.
1278 AddDiagonal(vertices,&newnumvertices,vindex,helpers[edgeIter->index],
1279 vertextypes, edgeTreeIterators, &edgeTree, helpers);
1280 }
1281 //helper(e j)�vi
1282 helpers[edgeIter->index] = vindex;
1283 }
1284 break;
1285 }
1286
1287 if(error) break;
1288 }
1289
1290 char *used = new char[newnumvertices];
1291 memset(used,0,newnumvertices*sizeof(char));
1292
1293 if(!error) {
1294 //return result
1295 long size;
1296 TPPLPoly mpoly;
1297 for(i=0;i<newnumvertices;i++) {
1298 if(used[i]) continue;
1299 v = &(vertices[i]);
1300 vnext = &(vertices[v->next]);
1301 size = 1;
1302 while(vnext!=v) {
1303 vnext = &(vertices[vnext->next]);
1304 size++;
1305 }
1306 mpoly.Init(size);
1307 v = &(vertices[i]);
1308 mpoly[0] = v->p;
1309 vnext = &(vertices[v->next]);
1310 size = 1;
1311 used[i] = 1;
1312 used[v->next] = 1;
1313 while(vnext!=v) {
1314 mpoly[size] = vnext->p;
1315 used[vnext->next] = 1;
1316 vnext = &(vertices[vnext->next]);
1317 size++;
1318 }
1319 monotonePolys->push_back(mpoly);
1320 }
1321 }
1322
1323 //cleanup
1324 delete [] vertices;
1325 delete [] priority;
1326 delete [] vertextypes;
1327 delete [] edgeTreeIterators;
1328 delete [] helpers;
1329 delete [] used;
1330
1331 if(error) {
1332 return 0;
1333 } else {
1334 return 1;
1335 }
1336 }
1337
1338 //adds a diagonal to the doubly-connected list of vertices
AddDiagonal(MonotoneVertex * vertices,long * numvertices,long index1,long index2,char * vertextypes,set<ScanLineEdge>::iterator * edgeTreeIterators,set<ScanLineEdge> * edgeTree,long * helpers)1339 void TPPLPartition::AddDiagonal(MonotoneVertex *vertices, long *numvertices, long index1, long index2,
1340 char *vertextypes, set<ScanLineEdge>::iterator *edgeTreeIterators,
1341 set<ScanLineEdge> *edgeTree, long *helpers)
1342 {
1343 long newindex1,newindex2;
1344
1345 newindex1 = *numvertices;
1346 (*numvertices)++;
1347 newindex2 = *numvertices;
1348 (*numvertices)++;
1349
1350 vertices[newindex1].p = vertices[index1].p;
1351 vertices[newindex2].p = vertices[index2].p;
1352
1353 vertices[newindex2].next = vertices[index2].next;
1354 vertices[newindex1].next = vertices[index1].next;
1355
1356 vertices[vertices[index2].next].previous = newindex2;
1357 vertices[vertices[index1].next].previous = newindex1;
1358
1359 vertices[index1].next = newindex2;
1360 vertices[newindex2].previous = index1;
1361
1362 vertices[index2].next = newindex1;
1363 vertices[newindex1].previous = index2;
1364
1365 //update all relevant structures
1366 vertextypes[newindex1] = vertextypes[index1];
1367 edgeTreeIterators[newindex1] = edgeTreeIterators[index1];
1368 helpers[newindex1] = helpers[index1];
1369 if(edgeTreeIterators[newindex1] != edgeTree->end())
1370 edgeTreeIterators[newindex1]->index = newindex1;
1371 vertextypes[newindex2] = vertextypes[index2];
1372 edgeTreeIterators[newindex2] = edgeTreeIterators[index2];
1373 helpers[newindex2] = helpers[index2];
1374 if(edgeTreeIterators[newindex2] != edgeTree->end())
1375 edgeTreeIterators[newindex2]->index = newindex2;
1376 }
1377
Below(TPPLPoint & p1,TPPLPoint & p2)1378 bool TPPLPartition::Below(TPPLPoint &p1, TPPLPoint &p2) {
1379 if(p1.y < p2.y) return true;
1380 else if(p1.y == p2.y) {
1381 if(p1.x < p2.x) return true;
1382 }
1383 return false;
1384 }
1385
1386 //sorts in the falling order of y values, if y is equal, x is used instead
operator ()(long index1,long index2)1387 bool TPPLPartition::VertexSorter::operator() (long index1, long index2) {
1388 if(vertices[index1].p.y > vertices[index2].p.y) return true;
1389 else if(vertices[index1].p.y == vertices[index2].p.y) {
1390 if(vertices[index1].p.x > vertices[index2].p.x) return true;
1391 }
1392 return false;
1393 }
1394
IsConvex(const TPPLPoint & p1,const TPPLPoint & p2,const TPPLPoint & p3) const1395 bool TPPLPartition::ScanLineEdge::IsConvex(const TPPLPoint& p1, const TPPLPoint& p2, const TPPLPoint& p3) const {
1396 tppl_float tmp;
1397 tmp = (p3.y-p1.y)*(p2.x-p1.x)-(p3.x-p1.x)*(p2.y-p1.y);
1398 if(tmp>0) return 1;
1399 else return 0;
1400 }
1401
operator <(const ScanLineEdge & other) const1402 bool TPPLPartition::ScanLineEdge::operator < (const ScanLineEdge & other) const {
1403 if(other.p1.y == other.p2.y) {
1404 if(p1.y == p2.y) {
1405 if(p1.y < other.p1.y) return true;
1406 else return false;
1407 }
1408 if(IsConvex(p1,p2,other.p1)) return true;
1409 else return false;
1410 } else if(p1.y == p2.y) {
1411 if(IsConvex(other.p1,other.p2,p1)) return false;
1412 else return true;
1413 } else if(p1.y < other.p1.y) {
1414 if(IsConvex(other.p1,other.p2,p1)) return false;
1415 else return true;
1416 } else {
1417 if(IsConvex(p1,p2,other.p1)) return true;
1418 else return false;
1419 }
1420 }
1421
1422 //triangulates monotone polygon
1423 //O(n) time, O(n) space complexity
TriangulateMonotone(TPPLPoly * inPoly,TPPLPolyList * triangles)1424 int TPPLPartition::TriangulateMonotone(TPPLPoly *inPoly, TPPLPolyList *triangles) {
1425 if(!inPoly->Valid()) return 0;
1426
1427 long i,i2,j,topindex,bottomindex,leftindex,rightindex,vindex;
1428 TPPLPoint *points = NULL;
1429 long numpoints;
1430 TPPLPoly triangle;
1431
1432 numpoints = inPoly->GetNumPoints();
1433 points = inPoly->GetPoints();
1434
1435 //trivial case
1436 if(numpoints == 3) {
1437 triangles->push_back(*inPoly);
1438 return 1;
1439 }
1440
1441 topindex = 0; bottomindex=0;
1442 for(i=1;i<numpoints;i++) {
1443 if(Below(points[i],points[bottomindex])) bottomindex = i;
1444 if(Below(points[topindex],points[i])) topindex = i;
1445 }
1446
1447 //check if the poly is really monotone
1448 i = topindex;
1449 while(i!=bottomindex) {
1450 i2 = i+1; if(i2>=numpoints) i2 = 0;
1451 if(!Below(points[i2],points[i])) return 0;
1452 i = i2;
1453 }
1454 i = bottomindex;
1455 while(i!=topindex) {
1456 i2 = i+1; if(i2>=numpoints) i2 = 0;
1457 if(!Below(points[i],points[i2])) return 0;
1458 i = i2;
1459 }
1460
1461 char *vertextypes = new char[numpoints];
1462 long *priority = new long[numpoints];
1463
1464 //merge left and right vertex chains
1465 priority[0] = topindex;
1466 vertextypes[topindex] = 0;
1467 leftindex = topindex+1; if(leftindex>=numpoints) leftindex = 0;
1468 rightindex = topindex-1; if(rightindex<0) rightindex = numpoints-1;
1469 for(i=1;i<(numpoints-1);i++) {
1470 if(leftindex==bottomindex) {
1471 priority[i] = rightindex;
1472 rightindex--; if(rightindex<0) rightindex = numpoints-1;
1473 vertextypes[priority[i]] = -1;
1474 } else if(rightindex==bottomindex) {
1475 priority[i] = leftindex;
1476 leftindex++; if(leftindex>=numpoints) leftindex = 0;
1477 vertextypes[priority[i]] = 1;
1478 } else {
1479 if(Below(points[leftindex],points[rightindex])) {
1480 priority[i] = rightindex;
1481 rightindex--; if(rightindex<0) rightindex = numpoints-1;
1482 vertextypes[priority[i]] = -1;
1483 } else {
1484 priority[i] = leftindex;
1485 leftindex++; if(leftindex>=numpoints) leftindex = 0;
1486 vertextypes[priority[i]] = 1;
1487 }
1488 }
1489 }
1490 priority[i] = bottomindex;
1491 vertextypes[bottomindex] = 0;
1492
1493 long *stack = new long[numpoints];
1494 long stackptr = 0;
1495
1496 stack[0] = priority[0];
1497 stack[1] = priority[1];
1498 stackptr = 2;
1499
1500 //for each vertex from top to bottom trim as many triangles as possible
1501 for(i=2;i<(numpoints-1);i++) {
1502 vindex = priority[i];
1503 if(vertextypes[vindex]!=vertextypes[stack[stackptr-1]]) {
1504 for(j=0;j<(stackptr-1);j++) {
1505 if(vertextypes[vindex]==1) {
1506 triangle.Triangle(points[stack[j+1]],points[stack[j]],points[vindex]);
1507 } else {
1508 triangle.Triangle(points[stack[j]],points[stack[j+1]],points[vindex]);
1509 }
1510 triangles->push_back(triangle);
1511 }
1512 stack[0] = priority[i-1];
1513 stack[1] = priority[i];
1514 stackptr = 2;
1515 } else {
1516 stackptr--;
1517 while(stackptr>0) {
1518 if(vertextypes[vindex]==1) {
1519 if(IsConvex(points[vindex],points[stack[stackptr-1]],points[stack[stackptr]])) {
1520 triangle.Triangle(points[vindex],points[stack[stackptr-1]],points[stack[stackptr]]);
1521 triangles->push_back(triangle);
1522 stackptr--;
1523 } else {
1524 break;
1525 }
1526 } else {
1527 if(IsConvex(points[vindex],points[stack[stackptr]],points[stack[stackptr-1]])) {
1528 triangle.Triangle(points[vindex],points[stack[stackptr]],points[stack[stackptr-1]]);
1529 triangles->push_back(triangle);
1530 stackptr--;
1531 } else {
1532 break;
1533 }
1534 }
1535 }
1536 stackptr++;
1537 stack[stackptr] = vindex;
1538 stackptr++;
1539 }
1540 }
1541 vindex = priority[i];
1542 for(j=0;j<(stackptr-1);j++) {
1543 if(vertextypes[stack[j+1]]==1) {
1544 triangle.Triangle(points[stack[j]],points[stack[j+1]],points[vindex]);
1545 } else {
1546 triangle.Triangle(points[stack[j+1]],points[stack[j]],points[vindex]);
1547 }
1548 triangles->push_back(triangle);
1549 }
1550
1551 delete [] priority;
1552 delete [] vertextypes;
1553 delete [] stack;
1554
1555 return 1;
1556 }
1557
Triangulate_MONO(TPPLPolyList * inpolys,TPPLPolyList * triangles)1558 int TPPLPartition::Triangulate_MONO(TPPLPolyList *inpolys, TPPLPolyList *triangles) {
1559 TPPLPolyList monotone;
1560 TPPLPolyList::iterator iter;
1561
1562 if(!MonotonePartition(inpolys,&monotone)) return 0;
1563 for(iter = monotone.begin(); iter!=monotone.end();iter++) {
1564 if(!TriangulateMonotone(&(*iter),triangles)) return 0;
1565 }
1566 return 1;
1567 }
1568
Triangulate_MONO(TPPLPoly * poly,TPPLPolyList * triangles)1569 int TPPLPartition::Triangulate_MONO(TPPLPoly *poly, TPPLPolyList *triangles) {
1570 TPPLPolyList polys;
1571 polys.push_back(*poly);
1572
1573 return Triangulate_MONO(&polys, triangles);
1574 }
1575