1 /*
2 * Poly2Tri Copyright (c) 2009-2018, Poly2Tri Contributors
3 * https://github.com/jhasse/poly2tri
4 *
5 * All rights reserved.
6 *
7 * Redistribution and use in source and binary forms, with or without modification,
8 * are permitted provided that the following conditions are met:
9 *
10 * * Redistributions of source code must retain the above copyright notice,
11 * this list of conditions and the following disclaimer.
12 * * Redistributions in binary form must reproduce the above copyright notice,
13 * this list of conditions and the following disclaimer in the documentation
14 * and/or other materials provided with the distribution.
15 * * Neither the name of Poly2Tri nor the names of its contributors may be
16 * used to endorse or promote products derived from this software without specific
17 * prior written permission.
18 *
19 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
20 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
21 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
22 * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
23 * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
24 * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
25 * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
26 * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
27 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
28 * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
29 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
30 */
31 #include "sweep.h"
32 #include "sweep_context.h"
33 #include "advancing_front.h"
34 #include "../common/utils.h"
35
36 #include <cassert>
37 #include <stdexcept>
38
39 namespace p2t {
40
41 // Triangulate simple polygon with holes
Triangulate(SweepContext & tcx)42 void Sweep::Triangulate(SweepContext& tcx)
43 {
44 tcx.InitTriangulation();
45 tcx.CreateAdvancingFront();
46 // Sweep points; build mesh
47 SweepPoints(tcx);
48 // Clean up
49 FinalizationPolygon(tcx);
50 }
51
SweepPoints(SweepContext & tcx)52 void Sweep::SweepPoints(SweepContext& tcx)
53 {
54 for (size_t i = 1; i < tcx.point_count(); i++) {
55 Point& point = *tcx.GetPoint(i);
56 Node* node = &PointEvent(tcx, point);
57 for (unsigned int j = 0; j < point.edge_list.size(); j++) {
58 EdgeEvent(tcx, point.edge_list[j], node);
59 }
60 }
61 }
62
FinalizationPolygon(SweepContext & tcx)63 void Sweep::FinalizationPolygon(SweepContext& tcx)
64 {
65 // Get an Internal triangle to start with
66 Triangle* t = tcx.front()->head()->next->triangle;
67 Point* p = tcx.front()->head()->next->point;
68 while (t && !t->GetConstrainedEdgeCW(*p)) {
69 t = t->NeighborCCW(*p);
70 }
71
72 // Collect interior triangles constrained by edges
73 if (t) {
74 tcx.MeshClean(*t);
75 }
76 }
77
PointEvent(SweepContext & tcx,Point & point)78 Node& Sweep::PointEvent(SweepContext& tcx, Point& point)
79 {
80 Node* node_ptr = tcx.LocateNode(point);
81 if (!node_ptr || !node_ptr->point || !node_ptr->next || !node_ptr->next->point)
82 {
83 throw std::runtime_error("PointEvent - null node");
84 }
85
86 Node& node = *node_ptr;
87 Node& new_node = NewFrontTriangle(tcx, point, node);
88
89 // Only need to check +epsilon since point never have smaller
90 // x value than node due to how we fetch nodes from the front
91 if (point.x <= node.point->x + EPSILON) {
92 Fill(tcx, node);
93 }
94
95 //tcx.AddNode(new_node);
96
97 FillAdvancingFront(tcx, new_node);
98 return new_node;
99 }
100
EdgeEvent(SweepContext & tcx,Edge * edge,Node * node)101 void Sweep::EdgeEvent(SweepContext& tcx, Edge* edge, Node* node)
102 {
103 tcx.edge_event.constrained_edge = edge;
104 tcx.edge_event.right = (edge->p->x > edge->q->x);
105
106 if (IsEdgeSideOfTriangle(*node->triangle, *edge->p, *edge->q)) {
107 return;
108 }
109
110 // For now we will do all needed filling
111 // TODO: integrate with flip process might give some better performance
112 // but for now this avoid the issue with cases that needs both flips and fills
113 FillEdgeEvent(tcx, edge, node);
114 EdgeEvent(tcx, *edge->p, *edge->q, node->triangle, *edge->q);
115 }
116
EdgeEvent(SweepContext & tcx,Point & ep,Point & eq,Triangle * triangle,Point & point)117 void Sweep::EdgeEvent(SweepContext& tcx, Point& ep, Point& eq, Triangle* triangle, Point& point)
118 {
119 if (triangle == nullptr) {
120 throw std::runtime_error("EdgeEvent - null triangle");
121 }
122 if (IsEdgeSideOfTriangle(*triangle, ep, eq)) {
123 return;
124 }
125
126 Point* p1 = triangle->PointCCW(point);
127 Orientation o1 = Orient2d(eq, *p1, ep);
128 if (o1 == COLLINEAR) {
129 if (triangle->Contains(&eq, p1)) {
130 triangle->MarkConstrainedEdge(&eq, p1);
131 // We are modifying the constraint maybe it would be better to
132 // not change the given constraint and just keep a variable for the new constraint
133 tcx.edge_event.constrained_edge->q = p1;
134 triangle = triangle->NeighborAcross(point);
135 EdgeEvent(tcx, ep, *p1, triangle, *p1);
136 } else {
137 throw std::runtime_error("EdgeEvent - collinear points not supported");
138 }
139 return;
140 }
141
142 Point* p2 = triangle->PointCW(point);
143 Orientation o2 = Orient2d(eq, *p2, ep);
144 if (o2 == COLLINEAR) {
145 if (triangle->Contains(&eq, p2)) {
146 triangle->MarkConstrainedEdge(&eq, p2);
147 // We are modifying the constraint maybe it would be better to
148 // not change the given constraint and just keep a variable for the new constraint
149 tcx.edge_event.constrained_edge->q = p2;
150 triangle = triangle->NeighborAcross(point);
151 EdgeEvent(tcx, ep, *p2, triangle, *p2);
152 } else {
153 throw std::runtime_error("EdgeEvent - collinear points not supported");
154 }
155 return;
156 }
157
158 if (o1 == o2) {
159 // Need to decide if we are rotating CW or CCW to get to a triangle
160 // that will cross edge
161 if (o1 == CW) {
162 triangle = triangle->NeighborCCW(point);
163 } else {
164 triangle = triangle->NeighborCW(point);
165 }
166 EdgeEvent(tcx, ep, eq, triangle, point);
167 } else {
168 // This triangle crosses constraint so lets flippin start!
169 assert(triangle);
170 FlipEdgeEvent(tcx, ep, eq, triangle, point);
171 }
172 }
173
IsEdgeSideOfTriangle(Triangle & triangle,Point & ep,Point & eq)174 bool Sweep::IsEdgeSideOfTriangle(Triangle& triangle, Point& ep, Point& eq)
175 {
176 const int index = triangle.EdgeIndex(&ep, &eq);
177
178 if (index != -1) {
179 triangle.MarkConstrainedEdge(index);
180 Triangle* t = triangle.GetNeighbor(index);
181 if (t) {
182 t->MarkConstrainedEdge(&ep, &eq);
183 }
184 return true;
185 }
186 return false;
187 }
188
NewFrontTriangle(SweepContext & tcx,Point & point,Node & node)189 Node& Sweep::NewFrontTriangle(SweepContext& tcx, Point& point, Node& node)
190 {
191 Triangle* triangle = new Triangle(point, *node.point, *node.next->point);
192
193 triangle->MarkNeighbor(*node.triangle);
194 tcx.AddToMap(triangle);
195
196 Node* new_node = new Node(point);
197 nodes_.push_back(new_node);
198
199 new_node->next = node.next;
200 new_node->prev = &node;
201 node.next->prev = new_node;
202 node.next = new_node;
203
204 if (!Legalize(tcx, *triangle)) {
205 tcx.MapTriangleToNodes(*triangle);
206 }
207
208 return *new_node;
209 }
210
Fill(SweepContext & tcx,Node & node)211 void Sweep::Fill(SweepContext& tcx, Node& node)
212 {
213 Triangle* triangle = new Triangle(*node.prev->point, *node.point, *node.next->point);
214
215 // TODO: should copy the constrained_edge value from neighbor triangles
216 // for now constrained_edge values are copied during the legalize
217 triangle->MarkNeighbor(*node.prev->triangle);
218 triangle->MarkNeighbor(*node.triangle);
219
220 tcx.AddToMap(triangle);
221
222 // Update the advancing front
223 node.prev->next = node.next;
224 node.next->prev = node.prev;
225
226 // If it was legalized the triangle has already been mapped
227 if (!Legalize(tcx, *triangle)) {
228 tcx.MapTriangleToNodes(*triangle);
229 }
230 }
231
FillAdvancingFront(SweepContext & tcx,Node & n)232 void Sweep::FillAdvancingFront(SweepContext& tcx, Node& n)
233 {
234
235 // Fill right holes
236 Node* node = n.next;
237
238 while (node && node->next) {
239 // if HoleAngle exceeds 90 degrees then break.
240 if (LargeHole_DontFill(node)) break;
241 Fill(tcx, *node);
242 node = node->next;
243 }
244
245 // Fill left holes
246 node = n.prev;
247
248 while (node && node->prev) {
249 // if HoleAngle exceeds 90 degrees then break.
250 if (LargeHole_DontFill(node)) break;
251 Fill(tcx, *node);
252 node = node->prev;
253 }
254
255 // Fill right basins
256 if (n.next && n.next->next) {
257 const double angle = BasinAngle(n);
258 if (angle < PI_3div4) {
259 FillBasin(tcx, n);
260 }
261 }
262 }
263
264 // True if HoleAngle exceeds 90 degrees.
LargeHole_DontFill(const Node * node) const265 bool Sweep::LargeHole_DontFill(const Node* node) const {
266
267 const Node* nextNode = node->next;
268 const Node* prevNode = node->prev;
269 if (!AngleExceeds90Degrees(node->point, nextNode->point, prevNode->point))
270 return false;
271
272 // Check additional points on front.
273 const Node* next2Node = nextNode->next;
274 // "..Plus.." because only want angles on same side as point being added.
275 if ((next2Node != nullptr) && !AngleExceedsPlus90DegreesOrIsNegative(node->point, next2Node->point, prevNode->point))
276 return false;
277
278 const Node* prev2Node = prevNode->prev;
279 // "..Plus.." because only want angles on same side as point being added.
280 if ((prev2Node != nullptr) && !AngleExceedsPlus90DegreesOrIsNegative(node->point, nextNode->point, prev2Node->point))
281 return false;
282
283 return true;
284 }
285
AngleExceeds90Degrees(const Point * origin,const Point * pa,const Point * pb) const286 bool Sweep::AngleExceeds90Degrees(const Point* origin, const Point* pa, const Point* pb) const {
287 const double angle = Angle(origin, pa, pb);
288 return ((angle > PI_div2) || (angle < -PI_div2));
289 }
290
AngleExceedsPlus90DegreesOrIsNegative(const Point * origin,const Point * pa,const Point * pb) const291 bool Sweep::AngleExceedsPlus90DegreesOrIsNegative(const Point* origin, const Point* pa, const Point* pb) const {
292 const double angle = Angle(origin, pa, pb);
293 return (angle > PI_div2) || (angle < 0);
294 }
295
Angle(const Point * origin,const Point * pa,const Point * pb) const296 double Sweep::Angle(const Point* origin, const Point* pa, const Point* pb) const {
297 /* Complex plane
298 * ab = cosA +i*sinA
299 * ab = (ax + ay*i)(bx + by*i) = (ax*bx + ay*by) + i(ax*by-ay*bx)
300 * atan2(y,x) computes the principal value of the argument function
301 * applied to the complex number x+iy
302 * Where x = ax*bx + ay*by
303 * y = ax*by - ay*bx
304 */
305 const double px = origin->x;
306 const double py = origin->y;
307 const double ax = pa->x - px;
308 const double ay = pa->y - py;
309 const double bx = pb->x - px;
310 const double by = pb->y - py;
311 const double x = ax * by - ay * bx;
312 const double y = ax * bx + ay * by;
313 return atan2(x, y);
314 }
315
BasinAngle(const Node & node) const316 double Sweep::BasinAngle(const Node& node) const
317 {
318 const double ax = node.point->x - node.next->next->point->x;
319 const double ay = node.point->y - node.next->next->point->y;
320 return atan2(ay, ax);
321 }
322
HoleAngle(const Node & node) const323 double Sweep::HoleAngle(const Node& node) const
324 {
325 /* Complex plane
326 * ab = cosA +i*sinA
327 * ab = (ax + ay*i)(bx + by*i) = (ax*bx + ay*by) + i(ax*by-ay*bx)
328 * atan2(y,x) computes the principal value of the argument function
329 * applied to the complex number x+iy
330 * Where x = ax*bx + ay*by
331 * y = ax*by - ay*bx
332 */
333 const double ax = node.next->point->x - node.point->x;
334 const double ay = node.next->point->y - node.point->y;
335 const double bx = node.prev->point->x - node.point->x;
336 const double by = node.prev->point->y - node.point->y;
337 return atan2(ax * by - ay * bx, ax * bx + ay * by);
338 }
339
Legalize(SweepContext & tcx,Triangle & t)340 bool Sweep::Legalize(SweepContext& tcx, Triangle& t)
341 {
342 // To legalize a triangle we start by finding if any of the three edges
343 // violate the Delaunay condition
344 for (int i = 0; i < 3; i++) {
345 if (t.delaunay_edge[i])
346 continue;
347
348 Triangle* ot = t.GetNeighbor(i);
349
350 if (ot) {
351 Point* p = t.GetPoint(i);
352 Point* op = ot->OppositePoint(t, *p);
353 int oi = ot->Index(op);
354
355 // If this is a Constrained Edge or a Delaunay Edge(only during recursive legalization)
356 // then we should not try to legalize
357 if (ot->constrained_edge[oi] || ot->delaunay_edge[oi]) {
358 t.constrained_edge[i] = ot->constrained_edge[oi];
359 continue;
360 }
361
362 bool inside = Incircle(*p, *t.PointCCW(*p), *t.PointCW(*p), *op);
363
364 if (inside) {
365 // Lets mark this shared edge as Delaunay
366 t.delaunay_edge[i] = true;
367 ot->delaunay_edge[oi] = true;
368
369 // Lets rotate shared edge one vertex CW to legalize it
370 RotateTrianglePair(t, *p, *ot, *op);
371
372 // We now got one valid Delaunay Edge shared by two triangles
373 // This gives us 4 new edges to check for Delaunay
374
375 // Make sure that triangle to node mapping is done only one time for a specific triangle
376 bool not_legalized = !Legalize(tcx, t);
377 if (not_legalized) {
378 tcx.MapTriangleToNodes(t);
379 }
380
381 not_legalized = !Legalize(tcx, *ot);
382 if (not_legalized)
383 tcx.MapTriangleToNodes(*ot);
384
385 // Reset the Delaunay edges, since they only are valid Delaunay edges
386 // until we add a new triangle or point.
387 // XXX: need to think about this. Can these edges be tried after we
388 // return to previous recursive level?
389 t.delaunay_edge[i] = false;
390 ot->delaunay_edge[oi] = false;
391
392 // If triangle have been legalized no need to check the other edges since
393 // the recursive legalization will handles those so we can end here.
394 return true;
395 }
396 }
397 }
398 return false;
399 }
400
Incircle(const Point & pa,const Point & pb,const Point & pc,const Point & pd) const401 bool Sweep::Incircle(const Point& pa, const Point& pb, const Point& pc, const Point& pd) const
402 {
403 const double adx = pa.x - pd.x;
404 const double ady = pa.y - pd.y;
405 const double bdx = pb.x - pd.x;
406 const double bdy = pb.y - pd.y;
407
408 const double adxbdy = adx * bdy;
409 const double bdxady = bdx * ady;
410 const double oabd = adxbdy - bdxady;
411
412 if (oabd <= 0)
413 return false;
414
415 const double cdx = pc.x - pd.x;
416 const double cdy = pc.y - pd.y;
417
418 const double cdxady = cdx * ady;
419 const double adxcdy = adx * cdy;
420 const double ocad = cdxady - adxcdy;
421
422 if (ocad <= 0)
423 return false;
424
425 const double bdxcdy = bdx * cdy;
426 const double cdxbdy = cdx * bdy;
427
428 const double alift = adx * adx + ady * ady;
429 const double blift = bdx * bdx + bdy * bdy;
430 const double clift = cdx * cdx + cdy * cdy;
431
432 const double det = alift * (bdxcdy - cdxbdy) + blift * ocad + clift * oabd;
433
434 return det > 0;
435 }
436
RotateTrianglePair(Triangle & t,Point & p,Triangle & ot,Point & op) const437 void Sweep::RotateTrianglePair(Triangle& t, Point& p, Triangle& ot, Point& op) const
438 {
439 Triangle* n1, *n2, *n3, *n4;
440 n1 = t.NeighborCCW(p);
441 n2 = t.NeighborCW(p);
442 n3 = ot.NeighborCCW(op);
443 n4 = ot.NeighborCW(op);
444
445 bool ce1, ce2, ce3, ce4;
446 ce1 = t.GetConstrainedEdgeCCW(p);
447 ce2 = t.GetConstrainedEdgeCW(p);
448 ce3 = ot.GetConstrainedEdgeCCW(op);
449 ce4 = ot.GetConstrainedEdgeCW(op);
450
451 bool de1, de2, de3, de4;
452 de1 = t.GetDelunayEdgeCCW(p);
453 de2 = t.GetDelunayEdgeCW(p);
454 de3 = ot.GetDelunayEdgeCCW(op);
455 de4 = ot.GetDelunayEdgeCW(op);
456
457 t.Legalize(p, op);
458 ot.Legalize(op, p);
459
460 // Remap delaunay_edge
461 ot.SetDelunayEdgeCCW(p, de1);
462 t.SetDelunayEdgeCW(p, de2);
463 t.SetDelunayEdgeCCW(op, de3);
464 ot.SetDelunayEdgeCW(op, de4);
465
466 // Remap constrained_edge
467 ot.SetConstrainedEdgeCCW(p, ce1);
468 t.SetConstrainedEdgeCW(p, ce2);
469 t.SetConstrainedEdgeCCW(op, ce3);
470 ot.SetConstrainedEdgeCW(op, ce4);
471
472 // Remap neighbors
473 // XXX: might optimize the markNeighbor by keeping track of
474 // what side should be assigned to what neighbor after the
475 // rotation. Now mark neighbor does lots of testing to find
476 // the right side.
477 t.ClearNeighbors();
478 ot.ClearNeighbors();
479 if (n1) ot.MarkNeighbor(*n1);
480 if (n2) t.MarkNeighbor(*n2);
481 if (n3) t.MarkNeighbor(*n3);
482 if (n4) ot.MarkNeighbor(*n4);
483 t.MarkNeighbor(ot);
484 }
485
FillBasin(SweepContext & tcx,Node & node)486 void Sweep::FillBasin(SweepContext& tcx, Node& node)
487 {
488 if (Orient2d(*node.point, *node.next->point, *node.next->next->point) == CCW) {
489 tcx.basin.left_node = node.next->next;
490 } else {
491 tcx.basin.left_node = node.next;
492 }
493
494 // Find the bottom and right node
495 tcx.basin.bottom_node = tcx.basin.left_node;
496 while (tcx.basin.bottom_node->next
497 && tcx.basin.bottom_node->point->y >= tcx.basin.bottom_node->next->point->y) {
498 tcx.basin.bottom_node = tcx.basin.bottom_node->next;
499 }
500 if (tcx.basin.bottom_node == tcx.basin.left_node) {
501 // No valid basin
502 return;
503 }
504
505 tcx.basin.right_node = tcx.basin.bottom_node;
506 while (tcx.basin.right_node->next
507 && tcx.basin.right_node->point->y < tcx.basin.right_node->next->point->y) {
508 tcx.basin.right_node = tcx.basin.right_node->next;
509 }
510 if (tcx.basin.right_node == tcx.basin.bottom_node) {
511 // No valid basins
512 return;
513 }
514
515 tcx.basin.width = tcx.basin.right_node->point->x - tcx.basin.left_node->point->x;
516 tcx.basin.left_highest = tcx.basin.left_node->point->y > tcx.basin.right_node->point->y;
517
518 FillBasinReq(tcx, tcx.basin.bottom_node);
519 }
520
FillBasinReq(SweepContext & tcx,Node * node)521 void Sweep::FillBasinReq(SweepContext& tcx, Node* node)
522 {
523 // if shallow stop filling
524 if (IsShallow(tcx, *node)) {
525 return;
526 }
527
528 Fill(tcx, *node);
529
530 if (node->prev == tcx.basin.left_node && node->next == tcx.basin.right_node) {
531 return;
532 } else if (node->prev == tcx.basin.left_node) {
533 Orientation o = Orient2d(*node->point, *node->next->point, *node->next->next->point);
534 if (o == CW) {
535 return;
536 }
537 node = node->next;
538 } else if (node->next == tcx.basin.right_node) {
539 Orientation o = Orient2d(*node->point, *node->prev->point, *node->prev->prev->point);
540 if (o == CCW) {
541 return;
542 }
543 node = node->prev;
544 } else {
545 // Continue with the neighbor node with lowest Y value
546 if (node->prev->point->y < node->next->point->y) {
547 node = node->prev;
548 } else {
549 node = node->next;
550 }
551 }
552
553 FillBasinReq(tcx, node);
554 }
555
IsShallow(SweepContext & tcx,Node & node)556 bool Sweep::IsShallow(SweepContext& tcx, Node& node)
557 {
558 double height;
559
560 if (tcx.basin.left_highest) {
561 height = tcx.basin.left_node->point->y - node.point->y;
562 } else {
563 height = tcx.basin.right_node->point->y - node.point->y;
564 }
565
566 // if shallow stop filling
567 if (tcx.basin.width > height) {
568 return true;
569 }
570 return false;
571 }
572
FillEdgeEvent(SweepContext & tcx,Edge * edge,Node * node)573 void Sweep::FillEdgeEvent(SweepContext& tcx, Edge* edge, Node* node)
574 {
575 if (tcx.edge_event.right) {
576 FillRightAboveEdgeEvent(tcx, edge, node);
577 } else {
578 FillLeftAboveEdgeEvent(tcx, edge, node);
579 }
580 }
581
FillRightAboveEdgeEvent(SweepContext & tcx,Edge * edge,Node * node)582 void Sweep::FillRightAboveEdgeEvent(SweepContext& tcx, Edge* edge, Node* node)
583 {
584 while (node->next->point->x < edge->p->x) {
585 // Check if next node is below the edge
586 if (Orient2d(*edge->q, *node->next->point, *edge->p) == CCW) {
587 FillRightBelowEdgeEvent(tcx, edge, *node);
588 } else {
589 node = node->next;
590 }
591 }
592 }
593
FillRightBelowEdgeEvent(SweepContext & tcx,Edge * edge,Node & node)594 void Sweep::FillRightBelowEdgeEvent(SweepContext& tcx, Edge* edge, Node& node)
595 {
596 if (node.point->x < edge->p->x) {
597 if (Orient2d(*node.point, *node.next->point, *node.next->next->point) == CCW) {
598 // Concave
599 FillRightConcaveEdgeEvent(tcx, edge, node);
600 } else {
601 // Convex
602 FillRightConvexEdgeEvent(tcx, edge, node);
603 // Retry this one
604 FillRightBelowEdgeEvent(tcx, edge, node);
605 }
606 }
607 }
608
FillRightConcaveEdgeEvent(SweepContext & tcx,Edge * edge,Node & node)609 void Sweep::FillRightConcaveEdgeEvent(SweepContext& tcx, Edge* edge, Node& node)
610 {
611 Fill(tcx, *node.next);
612 if (node.next->point != edge->p) {
613 // Next above or below edge?
614 if (Orient2d(*edge->q, *node.next->point, *edge->p) == CCW) {
615 // Below
616 if (Orient2d(*node.point, *node.next->point, *node.next->next->point) == CCW) {
617 // Next is concave
618 FillRightConcaveEdgeEvent(tcx, edge, node);
619 } else {
620 // Next is convex
621 }
622 }
623 }
624 }
625
FillRightConvexEdgeEvent(SweepContext & tcx,Edge * edge,Node & node)626 void Sweep::FillRightConvexEdgeEvent(SweepContext& tcx, Edge* edge, Node& node)
627 {
628 // Next concave or convex?
629 if (Orient2d(*node.next->point, *node.next->next->point, *node.next->next->next->point) == CCW) {
630 // Concave
631 FillRightConcaveEdgeEvent(tcx, edge, *node.next);
632 } else {
633 // Convex
634 // Next above or below edge?
635 if (Orient2d(*edge->q, *node.next->next->point, *edge->p) == CCW) {
636 // Below
637 FillRightConvexEdgeEvent(tcx, edge, *node.next);
638 } else {
639 // Above
640 }
641 }
642 }
643
FillLeftAboveEdgeEvent(SweepContext & tcx,Edge * edge,Node * node)644 void Sweep::FillLeftAboveEdgeEvent(SweepContext& tcx, Edge* edge, Node* node)
645 {
646 while (node->prev->point->x > edge->p->x) {
647 // Check if next node is below the edge
648 if (Orient2d(*edge->q, *node->prev->point, *edge->p) == CW) {
649 FillLeftBelowEdgeEvent(tcx, edge, *node);
650 } else {
651 node = node->prev;
652 }
653 }
654 }
655
FillLeftBelowEdgeEvent(SweepContext & tcx,Edge * edge,Node & node)656 void Sweep::FillLeftBelowEdgeEvent(SweepContext& tcx, Edge* edge, Node& node)
657 {
658 if (node.point->x > edge->p->x) {
659 if (Orient2d(*node.point, *node.prev->point, *node.prev->prev->point) == CW) {
660 // Concave
661 FillLeftConcaveEdgeEvent(tcx, edge, node);
662 } else {
663 // Convex
664 FillLeftConvexEdgeEvent(tcx, edge, node);
665 // Retry this one
666 FillLeftBelowEdgeEvent(tcx, edge, node);
667 }
668 }
669 }
670
FillLeftConvexEdgeEvent(SweepContext & tcx,Edge * edge,Node & node)671 void Sweep::FillLeftConvexEdgeEvent(SweepContext& tcx, Edge* edge, Node& node)
672 {
673 // Next concave or convex?
674 if (Orient2d(*node.prev->point, *node.prev->prev->point, *node.prev->prev->prev->point) == CW) {
675 // Concave
676 FillLeftConcaveEdgeEvent(tcx, edge, *node.prev);
677 } else {
678 // Convex
679 // Next above or below edge?
680 if (Orient2d(*edge->q, *node.prev->prev->point, *edge->p) == CW) {
681 // Below
682 FillLeftConvexEdgeEvent(tcx, edge, *node.prev);
683 } else {
684 // Above
685 }
686 }
687 }
688
FillLeftConcaveEdgeEvent(SweepContext & tcx,Edge * edge,Node & node)689 void Sweep::FillLeftConcaveEdgeEvent(SweepContext& tcx, Edge* edge, Node& node)
690 {
691 Fill(tcx, *node.prev);
692 if (node.prev->point != edge->p) {
693 // Next above or below edge?
694 if (Orient2d(*edge->q, *node.prev->point, *edge->p) == CW) {
695 // Below
696 if (Orient2d(*node.point, *node.prev->point, *node.prev->prev->point) == CW) {
697 // Next is concave
698 FillLeftConcaveEdgeEvent(tcx, edge, node);
699 } else {
700 // Next is convex
701 }
702 }
703 }
704 }
705
FlipEdgeEvent(SweepContext & tcx,Point & ep,Point & eq,Triangle * t,Point & p)706 void Sweep::FlipEdgeEvent(SweepContext& tcx, Point& ep, Point& eq, Triangle* t, Point& p)
707 {
708 assert(t);
709 Triangle* ot_ptr = t->NeighborAcross(p);
710 if (ot_ptr == nullptr)
711 {
712 throw std::runtime_error("FlipEdgeEvent - null neighbor across");
713 }
714 Triangle& ot = *ot_ptr;
715 Point& op = *ot.OppositePoint(*t, p);
716
717 if (InScanArea(p, *t->PointCCW(p), *t->PointCW(p), op)) {
718 // Lets rotate shared edge one vertex CW
719 RotateTrianglePair(*t, p, ot, op);
720 tcx.MapTriangleToNodes(*t);
721 tcx.MapTriangleToNodes(ot);
722
723 if (p == eq && op == ep) {
724 if (eq == *tcx.edge_event.constrained_edge->q && ep == *tcx.edge_event.constrained_edge->p) {
725 t->MarkConstrainedEdge(&ep, &eq);
726 ot.MarkConstrainedEdge(&ep, &eq);
727 Legalize(tcx, *t);
728 Legalize(tcx, ot);
729 } else {
730 // XXX: I think one of the triangles should be legalized here?
731 }
732 } else {
733 Orientation o = Orient2d(eq, op, ep);
734 t = &NextFlipTriangle(tcx, (int)o, *t, ot, p, op);
735 FlipEdgeEvent(tcx, ep, eq, t, p);
736 }
737 } else {
738 Point& newP = NextFlipPoint(ep, eq, ot, op);
739 FlipScanEdgeEvent(tcx, ep, eq, *t, ot, newP);
740 EdgeEvent(tcx, ep, eq, t, p);
741 }
742 }
743
NextFlipTriangle(SweepContext & tcx,int o,Triangle & t,Triangle & ot,Point & p,Point & op)744 Triangle& Sweep::NextFlipTriangle(SweepContext& tcx, int o, Triangle& t, Triangle& ot, Point& p, Point& op)
745 {
746 if (o == CCW) {
747 // ot is not crossing edge after flip
748 int edge_index = ot.EdgeIndex(&p, &op);
749 ot.delaunay_edge[edge_index] = true;
750 Legalize(tcx, ot);
751 ot.ClearDelunayEdges();
752 return t;
753 }
754
755 // t is not crossing edge after flip
756 int edge_index = t.EdgeIndex(&p, &op);
757
758 t.delaunay_edge[edge_index] = true;
759 Legalize(tcx, t);
760 t.ClearDelunayEdges();
761 return ot;
762 }
763
NextFlipPoint(Point & ep,Point & eq,Triangle & ot,Point & op)764 Point& Sweep::NextFlipPoint(Point& ep, Point& eq, Triangle& ot, Point& op)
765 {
766 Orientation o2d = Orient2d(eq, op, ep);
767 if (o2d == CW) {
768 // Right
769 return *ot.PointCCW(op);
770 } else if (o2d == CCW) {
771 // Left
772 return *ot.PointCW(op);
773 }
774 throw std::runtime_error("[Unsupported] Opposing point on constrained edge");
775 }
776
FlipScanEdgeEvent(SweepContext & tcx,Point & ep,Point & eq,Triangle & flip_triangle,Triangle & t,Point & p)777 void Sweep::FlipScanEdgeEvent(SweepContext& tcx, Point& ep, Point& eq, Triangle& flip_triangle,
778 Triangle& t, Point& p)
779 {
780 Triangle* ot_ptr = t.NeighborAcross(p);
781 if (ot_ptr == nullptr) {
782 throw std::runtime_error("FlipScanEdgeEvent - null neighbor across");
783 }
784
785 Point* op_ptr = ot_ptr->OppositePoint(t, p);
786 if (op_ptr == nullptr) {
787 throw std::runtime_error("FlipScanEdgeEvent - null opposing point");
788 }
789
790 Point* p1 = flip_triangle.PointCCW(eq);
791 Point* p2 = flip_triangle.PointCW(eq);
792 if (p1 == nullptr || p2 == nullptr) {
793 throw std::runtime_error("FlipScanEdgeEvent - null on either of points");
794 }
795
796 Triangle& ot = *ot_ptr;
797 Point& op = *op_ptr;
798
799 if (InScanArea(eq, *p1, *p2, op)) {
800 // flip with new edge op->eq
801 FlipEdgeEvent(tcx, eq, op, &ot, op);
802 // TODO: Actually I just figured out that it should be possible to
803 // improve this by getting the next ot and op before the the above
804 // flip and continue the flipScanEdgeEvent here
805 // set new ot and op here and loop back to inScanArea test
806 // also need to set a new flip_triangle first
807 // Turns out at first glance that this is somewhat complicated
808 // so it will have to wait.
809 } else {
810 Point& newP = NextFlipPoint(ep, eq, ot, op);
811 FlipScanEdgeEvent(tcx, ep, eq, flip_triangle, ot, newP);
812 }
813 }
814
~Sweep()815 Sweep::~Sweep() {
816
817 // Clean up memory
818 for(size_t i = 0; i < nodes_.size(); i++) {
819 delete nodes_[i];
820 }
821
822 }
823
824 } // namespace p2t
825