1 2 /** D port (betterC) of the earcut polygon triangulation library. 3 4 Ported from: https://github.com/mapbox/earcut.hpp 5 6 Copyright: 7 Copyright (c) 2020, Ferhat Kurtulmuş. 8 License: 9 $(HTTP www.boost.org/LICENSE_1_0.txt, Boost License 1.0). 10 */ 11 12 module earcutd; 13 14 version(LDC){ 15 version(D_BetterC){ 16 pragma(LDC_no_moduleinfo); 17 } 18 } 19 20 import dvector; 21 22 private { 23 import core.stdc.stdlib; 24 import std.traits; 25 import std.math; 26 import std.stdint; 27 import std.algorithm.comparison: min, max; 28 } 29 30 struct Earcut(N, Polygon) { 31 32 alias ASeq = TemplateArgsOf!Polygon; 33 alias VecPoint = ASeq[0]; 34 alias ASeq2 = TemplateArgsOf!VecPoint; 35 alias Point = ASeq2[0]; 36 37 Dvector!N indices; 38 size_t vertices = 0; 39 40 struct Node { 41 42 N i; 43 double x; 44 double y; 45 46 @nogc nothrow: 47 int opCmp(Node rhs) { 48 if (x < rhs.x) return -1; 49 if (rhs.x < x) return 1; 50 return 0; 51 } 52 53 // previous and next vertice nodes in a polygon ring 54 Node* prev; 55 Node* next; 56 57 // z-order curve value 58 int32_t z; 59 60 // previous and next nodes in z-order 61 Node* prevZ; 62 Node* nextZ; 63 64 // indicates whether this is a steiner point 65 bool steiner = false; 66 67 this(N index, double x_, double y_){ 68 i = index; 69 x = x_; 70 y = y_; 71 } 72 } 73 74 bool hashing; 75 double minX, maxX; 76 double minY, maxY; 77 double inv_size = 0.0; 78 79 struct ObjectPool(T) { 80 @nogc nothrow: 81 public: 82 this(size_t blockSize_) { 83 reset(blockSize_); 84 } 85 ~this() { 86 clear(); 87 } 88 89 T* construct(N)(N i, double x, double y) { 90 if (currentIndex >= blockSize) { 91 currentBlock = cast(T*)malloc(blockSize*T.sizeof); 92 allocations.pushBack(currentBlock); 93 currentIndex = 0; 94 } 95 T* object = ¤tBlock[currentIndex++]; 96 *object = T(i, x, y); 97 return object; 98 } 99 void reset(size_t newBlockSize) { 100 foreach (allocation; allocations) { 101 core.stdc.stdlib.free(allocation); 102 } 103 allocations.free(); 104 blockSize = max(1, newBlockSize); 105 currentBlock = null; 106 currentIndex = blockSize; 107 } 108 void clear() { reset(blockSize); } 109 private: 110 T* currentBlock = null; 111 size_t currentIndex = 1; 112 size_t blockSize = 1; 113 Dvector!(T*) allocations; 114 } 115 116 ObjectPool!Node nodes; 117 118 @nogc nothrow: 119 120 void run(ref Polygon points){ 121 // reset 122 indices.free; 123 vertices = 0; 124 125 if (!points.length) return; 126 127 double x; 128 double y; 129 int threshold = 80; 130 size_t len = 0; 131 132 for (size_t i = 0; threshold >= 0 && i < points.length; i++) { 133 threshold -= cast(int)(points[i].length); 134 len += points[i].length; 135 } 136 137 //estimate size of nodes and indices 138 nodes.reset(len * 3 / 2); 139 indices.reserve(len + points[0].length); 140 141 Node* outerNode = linkedList(points[0], true); 142 if (!outerNode || outerNode.prev == outerNode.next) return; 143 144 if (points.length > 1) outerNode = eliminateHoles(points, outerNode); 145 146 // if the shape is not too simple, we'll use z-order curve hash later; calculate polygon bbox 147 hashing = threshold < 0; 148 if (hashing) { 149 Node* p = outerNode.next; 150 minX = maxX = outerNode.x; 151 minY = maxY = outerNode.y; 152 do { 153 x = p.x; 154 y = p.y; 155 minX = min(minX, x); 156 minY = min(minY, y); 157 maxX = max(maxX, x); 158 maxY = max(maxY, y); 159 p = p.next; 160 } while (p != outerNode); 161 162 // minX, minY and size are later used to transform coords into integers for z-order calculation 163 inv_size = max(maxX - minX, maxY - minY); 164 inv_size = inv_size != .0 ? (1. / inv_size) : .0; 165 } 166 167 earcutLinked(outerNode); 168 169 nodes.clear(); 170 } 171 172 Node* linkedList(Ring)(ref Ring points, const bool clockwise) { 173 //using Point = typename Ring::value_type; 174 double sum = 0; 175 const size_t len = points.length; 176 size_t i, j; 177 Node* last = null; 178 179 // calculate original winding order of a polygon ring 180 for (i = 0, j = len > 0 ? len - 1 : 0; i < len; j = i++) { 181 const p1 = points[i]; 182 const p2 = points[j]; 183 const double p20 = p2[0]; 184 const double p10 = p1[0]; 185 const double p11 = p1[1]; 186 const double p21 = p2[1]; 187 188 sum += (p20 - p10) * (p11 + p21); 189 } 190 191 // link points into circular doubly-linked list in the specified winding order 192 if (clockwise == (sum > 0)) { 193 for (i = 0; i < len; i++) last = insertNode(vertices + i, points[i], last); 194 } else { 195 for (i = len; i-- > 0;) last = insertNode(vertices + i, points[i], last); 196 } 197 198 if (last && equals(last, last.next)) { 199 removeNode(last); 200 last = last.next; 201 } 202 203 vertices += len; 204 205 return last; 206 } 207 208 Node* filterPoints(Node* start, Node* end = null){ 209 if (!end) end = start; 210 211 Node* p = start; 212 bool again; 213 do { 214 again = false; 215 216 if (!p.steiner && (equals(p, p.next) || area(p.prev, p, p.next) == 0)) { 217 removeNode(p); 218 p = end = p.prev; 219 220 if (p == p.next) break; 221 again = true; 222 223 } else { 224 p = p.next; 225 } 226 } while (again || p != end); 227 228 return end; 229 } 230 231 void earcutLinked(Node* ear, int pass = 0){ 232 if (!ear) return; 233 234 // interlink polygon nodes in z-order 235 if (!pass && hashing) indexCurve(ear); 236 237 Node* stop = ear; 238 Node* prev; 239 Node* next; 240 241 int iterations = 0; 242 243 // iterate through ears, slicing them one by one 244 while (ear.prev != ear.next) { 245 iterations++; 246 prev = ear.prev; 247 next = ear.next; 248 249 if (hashing ? isEarHashed(ear) : isEar(ear)) { 250 // cut off the triangle 251 indices.pushBack(prev.i); 252 indices.pushBack(ear.i); 253 indices.pushBack(next.i); 254 255 removeNode(ear); 256 257 // skipping the next vertice leads to less sliver triangles 258 ear = next.next; 259 stop = next.next; 260 261 continue; 262 } 263 264 ear = next; 265 266 // if we looped through the whole remaining polygon and can't find any more ears 267 if (ear == stop) { 268 // try filtering points and slicing again 269 if (!pass) earcutLinked(filterPoints(ear), 1); 270 271 // if this didn't work, try curing all small self-intersections locally 272 else if (pass == 1) { 273 ear = cureLocalIntersections(filterPoints(ear)); 274 earcutLinked(ear, 2); 275 276 // as a last resort, try splitting the remaining polygon into two 277 } else if (pass == 2) splitEarcut(ear); 278 279 break; 280 } 281 } 282 283 } 284 285 bool isEar(Node* ear) { 286 Node* a = ear.prev; 287 Node* b = ear; 288 Node* c = ear.next; 289 290 if (area(a, b, c) >= 0) return false; // reflex, can't be an ear 291 292 // now make sure we don't have other points inside the potential ear 293 Node* p = ear.next.next; 294 295 while (p != ear.prev) { 296 if (pointInTriangle(a.x, a.y, b.x, b.y, c.x, c.y, p.x, p.y) && 297 area(p.prev, p, p.next) >= 0) return false; 298 p = p.next; 299 } 300 301 return true; 302 } 303 304 bool isEarHashed(Node* ear) { 305 Node* a = ear.prev; 306 Node* b = ear; 307 Node* c = ear.next; 308 309 if (area(a, b, c) >= 0) return false; // reflex, can't be an ear 310 311 // triangle bbox; min & max are calculated like this for speed 312 const double minTX = min(a.x, min(b.x, c.x)); 313 const double minTY = min(a.y, min(b.y, c.y)); 314 const double maxTX = max(a.x, max(b.x, c.x)); 315 const double maxTY = max(a.y, max(b.y, c.y)); 316 317 // z-order range for the current triangle bbox; 318 const int32_t minZ = zOrder(minTX, minTY); 319 const int32_t maxZ = zOrder(maxTX, maxTY); 320 321 // first look for points inside the triangle in increasing z-order 322 Node* p = ear.nextZ; 323 324 while (p && p.z <= maxZ) { 325 if (p != ear.prev && p != ear.next && 326 pointInTriangle(a.x, a.y, b.x, b.y, c.x, c.y, p.x, p.y) && 327 area(p.prev, p, p.next) >= 0) return false; 328 p = p.nextZ; 329 } 330 331 // then look for points in decreasing z-order 332 p = ear.prevZ; 333 334 while (p && p.z >= minZ) { 335 if (p != ear.prev && p != ear.next && 336 pointInTriangle(a.x, a.y, b.x, b.y, c.x, c.y, p.x, p.y) && 337 area(p.prev, p, p.next) >= 0) return false; 338 p = p.prevZ; 339 } 340 341 return true; 342 } 343 344 Node* cureLocalIntersections(Node* start) { 345 Node* p = start; 346 do { 347 Node* a = p.prev; 348 Node* b = p.next.next; 349 350 // a self-intersection where edge (v[i-1],v[i]) intersects (v[i+1],v[i+2]) 351 if (!equals(a, b) && intersects(a, p, p.next, b) && locallyInside(a, b) && locallyInside(b, a)) { 352 indices.pushBack(a.i); 353 indices.pushBack(p.i); 354 indices.pushBack(b.i); 355 356 // remove two nodes involved 357 removeNode(p); 358 removeNode(p.next); 359 360 p = start = b; 361 } 362 p = p.next; 363 } while (p != start); 364 365 return filterPoints(p); 366 } 367 368 void splitEarcut(Node* start) { 369 // look for a valid diagonal that divides the polygon into two 370 Node* a = start; 371 do { 372 Node* b = a.next.next; 373 while (b != a.prev) { 374 if (a.i != b.i && isValidDiagonal(a, b)) { 375 // split the polygon in two by the diagonal 376 Node* c = splitPolygon(a, b); 377 378 // filter colinear points around the cuts 379 a = filterPoints(a, a.next); 380 c = filterPoints(c, c.next); 381 382 // run earcut on each half 383 earcutLinked(a); 384 earcutLinked(c); 385 return; 386 } 387 b = b.next; 388 } 389 a = a.next; 390 } while (a != start); 391 } 392 393 Node* eliminateHoles(Polygon)(ref Polygon points, Node* outerNode){ 394 const size_t len = points.length; 395 396 Dvector!(Node*) queue; 397 for (size_t i = 1; i < len; i++) { 398 Node* list = linkedList(points[i], false); 399 if (list) { 400 if (list == list.next) list.steiner = true; 401 queue.pushBack(getLeftmost(list)); 402 } 403 } 404 405 import std.algorithm.sorting; 406 407 auto _q = queue.slice.sort!("a < b"); 408 409 // process holes from left to right 410 for (size_t i = 0; i < queue.length; i++) { 411 eliminateHole(_q[i], outerNode); 412 outerNode = filterPoints(outerNode, outerNode.next); 413 } 414 415 queue.free; 416 return outerNode; 417 } 418 419 void eliminateHole(Node* hole, Node* outerNode) { 420 outerNode = findHoleBridge(hole, outerNode); 421 if (outerNode) { 422 Node* b = splitPolygon(outerNode, hole); 423 424 // filter out colinear points around cuts 425 filterPoints(outerNode, outerNode.next); 426 filterPoints(b, b.next); 427 } 428 } 429 430 // David Eberly's algorithm for finding a bridge between hole and outer polygon 431 Node* findHoleBridge(Node* hole, Node* outerNode) { 432 Node* p = outerNode; 433 double hx = hole.x; 434 double hy = hole.y; 435 double qx = -double.max; 436 Node* m = null; 437 438 // find a segment intersected by a ray from the hole's leftmost Vertex to the left; 439 // segment's endpoint with lesser x will be potential connection Vertex 440 do { 441 if (hy <= p.y && hy >= p.next.y && p.next.y != p.y) { 442 double x = p.x + (hy - p.y) * (p.next.x - p.x) / (p.next.y - p.y); 443 if (x <= hx && x > qx) { 444 qx = x; 445 if (x == hx) { 446 if (hy == p.y) return p; 447 if (hy == p.next.y) return p.next; 448 } 449 m = p.x < p.next.x ? p : p.next; 450 } 451 } 452 p = p.next; 453 } while (p != outerNode); 454 455 if (!m) return null; 456 457 if (hx == qx) return m; // hole touches outer segment; pick leftmost endpoint 458 459 // look for points inside the triangle of hole Vertex, segment intersection and endpoint; 460 // if there are no points found, we have a valid connection; 461 // otherwise choose the Vertex of the minimum angle with the ray as connection Vertex 462 463 Node* stop = m; 464 double tanMin = double.max; 465 double tanCur = 0; 466 467 p = m; 468 double mx = m.x; 469 double my = m.y; 470 471 do { 472 if (hx >= p.x && p.x >= mx && hx != p.x && 473 pointInTriangle(hy < my ? hx : qx, hy, mx, my, hy < my ? qx : hx, hy, p.x, p.y)) { 474 475 tanCur = fabs(hy - p.y) / (hx - p.x); // tangential 476 477 if (locallyInside(p, hole) && 478 (tanCur < tanMin || (tanCur == tanMin && (p.x > m.x || sectorContainsSector(m, p))))) { 479 m = p; 480 tanMin = tanCur; 481 } 482 } 483 484 p = p.next; 485 } while (p != stop); 486 487 return m; 488 } 489 490 bool sectorContainsSector(Node* m, Node* p) { 491 return area(m.prev, m, p.prev) < 0 && area(p.next, m, m.next) < 0; 492 } 493 494 void indexCurve(Node* start) { 495 assert(start); 496 Node* p = start; 497 498 do { 499 p.z = p.z ? p.z : zOrder(p.x, p.y); 500 p.prevZ = p.prev; 501 p.nextZ = p.next; 502 p = p.next; 503 } while (p != start); 504 505 p.prevZ.nextZ = null; 506 p.prevZ = null; 507 508 sortLinked(p); 509 } 510 511 Node* sortLinked(Node* list) { 512 assert(list); 513 Node* p; 514 Node* q; 515 Node* e; 516 Node* tail; 517 int i, numMerges, pSize, qSize; 518 int inSize = 1; 519 520 for (;;) { 521 p = list; 522 list = null; 523 tail = null; 524 numMerges = 0; 525 526 while (p) { 527 numMerges++; 528 q = p; 529 pSize = 0; 530 for (i = 0; i < inSize; i++) { 531 pSize++; 532 q = q.nextZ; 533 if (!q) break; 534 } 535 536 qSize = inSize; 537 538 while (pSize > 0 || (qSize > 0 && q)) { 539 540 if (pSize == 0) { 541 e = q; 542 q = q.nextZ; 543 qSize--; 544 } else if (qSize == 0 || !q) { 545 e = p; 546 p = p.nextZ; 547 pSize--; 548 } else if (p.z <= q.z) { 549 e = p; 550 p = p.nextZ; 551 pSize--; 552 } else { 553 e = q; 554 q = q.nextZ; 555 qSize--; 556 } 557 558 if (tail) tail.nextZ = e; 559 else list = e; 560 561 e.prevZ = tail; 562 tail = e; 563 } 564 565 p = q; 566 } 567 568 tail.nextZ = null; 569 570 if (numMerges <= 1) return list; 571 572 inSize *= 2; 573 } 574 } 575 576 int32_t zOrder(const double x_, const double y_) { 577 // coords are transformed into non-negative 15-bit integer range 578 int32_t x = cast(int32_t)(32767.0 * (x_ - minX) * inv_size); 579 int32_t y = cast(int32_t)(32767.0 * (y_ - minY) * inv_size); 580 581 x = (x | (x << 8)) & 0x00FF00FF; 582 x = (x | (x << 4)) & 0x0F0F0F0F; 583 x = (x | (x << 2)) & 0x33333333; 584 x = (x | (x << 1)) & 0x55555555; 585 586 y = (y | (y << 8)) & 0x00FF00FF; 587 y = (y | (y << 4)) & 0x0F0F0F0F; 588 y = (y | (y << 2)) & 0x33333333; 589 y = (y | (y << 1)) & 0x55555555; 590 591 return x | (y << 1); 592 } 593 594 Node* getLeftmost(Node* start) { 595 Node* p = start; 596 Node* leftmost = start; 597 do { 598 if (p.x < leftmost.x || (p.x == leftmost.x && p.y < leftmost.y)) 599 leftmost = p; 600 p = p.next; 601 } while (p != start); 602 603 return leftmost; 604 } 605 606 bool pointInTriangle(double ax, double ay, double bx, double by, double cx, double cy, double px, double py) const { 607 return (cx - px) * (ay - py) - (ax - px) * (cy - py) >= 0 && 608 (ax - px) * (by - py) - (bx - px) * (ay - py) >= 0 && 609 (bx - px) * (cy - py) - (cx - px) * (by - py) >= 0; 610 } 611 612 bool isValidDiagonal(Node* a, Node* b) { 613 return a.next.i != b.i && a.prev.i != b.i && !intersectsPolygon(a, b) && // dones't intersect other edges 614 ((locallyInside(a, b) && locallyInside(b, a) && middleInside(a, b) && // locally visible 615 (area(a.prev, a, b.prev) != 0.0 || area(a, b.prev, b) != 0.0)) || // does not create opposite-facing sectors 616 (equals(a, b) && area(a.prev, a, a.next) > 0 && area(b.prev, b, b.next) > 0)); // special zero-length case 617 } 618 619 double area(Node* p, Node* q, Node* r) const { 620 return (q.y - p.y) * (r.x - q.x) - (q.x - p.x) * (r.y - q.y); 621 } 622 623 bool equals(Node* p1, Node* p2) { 624 return p1.x == p2.x && p1.y == p2.y; 625 } 626 627 bool intersects(Node* p1, Node* q1, Node* p2, Node* q2) { 628 int o1 = sign(area(p1, q1, p2)); 629 int o2 = sign(area(p1, q1, q2)); 630 int o3 = sign(area(p2, q2, p1)); 631 int o4 = sign(area(p2, q2, q1)); 632 633 if (o1 != o2 && o3 != o4) return true; // general case 634 635 if (o1 == 0 && onSegment(p1, p2, q1)) return true; // p1, q1 and p2 are collinear and p2 lies on p1q1 636 if (o2 == 0 && onSegment(p1, q2, q1)) return true; // p1, q1 and q2 are collinear and q2 lies on p1q1 637 if (o3 == 0 && onSegment(p2, p1, q2)) return true; // p2, q2 and p1 are collinear and p1 lies on p2q2 638 if (o4 == 0 && onSegment(p2, q1, q2)) return true; // p2, q2 and q1 are collinear and q1 lies on p2q2 639 640 return false; 641 } 642 643 bool onSegment(Node* p, Node* q, Node* r) { 644 return q.x <= max(p.x, r.x) && 645 q.x >= min(p.x, r.x) && 646 q.y <= max(p.y, r.y) && 647 q.y >= min(p.y, r.y); 648 } 649 650 int sign(double val) { 651 return (0.0 < val) - (val < 0.0); 652 } 653 654 bool intersectsPolygon(Node* a, Node* b) { 655 Node* p = a; 656 do { 657 if (p.i != a.i && p.next.i != a.i && p.i != b.i && p.next.i != b.i && 658 intersects(p, p.next, a, b)) return true; 659 p = p.next; 660 } while (p != a); 661 662 return false; 663 } 664 665 bool locallyInside(Node* a, Node* b) { 666 return area(a.prev, a, a.next) < 0 ? 667 area(a, b, a.next) >= 0 && area(a, a.prev, b) >= 0 : 668 area(a, b, a.prev) < 0 || area(a, a.next, b) < 0; 669 } 670 671 bool middleInside(Node* a, Node* b) { 672 Node* p = a; 673 bool inside = false; 674 double px = (a.x + b.x) / 2; 675 double py = (a.y + b.y) / 2; 676 do { 677 if (((p.y > py) != (p.next.y > py)) && p.next.y != p.y && 678 (px < (p.next.x - p.x) * (py - p.y) / (p.next.y - p.y) + p.x)) 679 inside = !inside; 680 p = p.next; 681 } while (p != a); 682 683 return inside; 684 } 685 686 Node* splitPolygon(Node* a, Node* b) { 687 Node* a2 = nodes.construct(a.i, a.x, a.y); 688 Node* b2 = nodes.construct(b.i, b.x, b.y); 689 Node* an = a.next; 690 Node* bp = b.prev; 691 692 a.next = b; 693 b.prev = a; 694 695 a2.next = an; 696 an.prev = a2; 697 698 b2.next = a2; 699 a2.prev = b2; 700 701 bp.next = b2; 702 b2.prev = bp; 703 704 return b2; 705 } 706 707 Node* insertNode(size_t i, ref Point pt, Node* last) { 708 Node* p = nodes.construct(cast(N)i, pt[0], pt[1]); 709 710 if (!last) { 711 p.prev = p; 712 p.next = p; 713 714 } else { 715 assert(last); 716 p.next = last.next; 717 p.prev = last; 718 last.next.prev = p; 719 last.next = p; 720 } 721 return p; 722 } 723 724 void removeNode(Node* p) { 725 p.next.prev = p.prev; 726 p.prev.next = p.next; 727 728 if (p.prevZ) p.prevZ.nextZ = p.nextZ; 729 if (p.nextZ) p.nextZ.prevZ = p.prevZ; 730 } 731 732 } 733 734 // those are optional. You can just use Tuple!(int, int) 735 struct Pair(T){ 736 T x; 737 T y; 738 @nogc nothrow: 739 inout(T) opIndex(size_t index) inout { 740 T[2] tmp = [x, y]; 741 return tmp[index]; 742 } 743 744 void opIndexAssign(T)(T value, size_t index){ 745 T*[2] tmp = [&x, &y]; 746 *tmp[i] = value; 747 } 748 } 749 750 struct Pair2(T){ 751 T[2] coord; 752 753 alias coord this; 754 755 this(T x, T y,){ 756 coord[0] = x; 757 coord[1] = y; 758 } 759 }