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