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