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