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