/[zanavi_public1]/navit/navit/maptool/p2t/sweep/sweep.c
ZANavi

Contents of /navit/navit/maptool/p2t/sweep/sweep.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 37 - (show annotations) (download)
Sat Mar 8 21:37:20 2014 UTC (10 years ago) by zoff99
File MIME type: text/plain
File size: 31927 byte(s)
new market version, lots of new features
1 /*
2 * This file is a part of the C port of the Poly2Tri library
3 * Porting to C done by (c) Barak Itkin <lightningismyname@gmail.com>
4 * http://code.google.com/p/poly2tri-c/
5 *
6 * Poly2Tri Copyright (c) 2009-2010, Poly2Tri Contributors
7 * http://code.google.com/p/poly2tri/
8 *
9 * All rights reserved.
10 *
11 * Redistribution and use in source and binary forms, with or without modification,
12 * are permitted provided that the following conditions are met:
13 *
14 * * Redistributions of source code must retain the above copyright notice,
15 * this list of conditions and the following disclaimer.
16 * * Redistributions in binary form must reproduce the above copyright notice,
17 * this list of conditions and the following disclaimer in the documentation
18 * and/or other materials provided with the distribution.
19 * * Neither the name of Poly2Tri nor the names of its contributors may be
20 * used to endorse or promote products derived from this software without specific
21 * prior written permission.
22 *
23 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
24 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
25 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
26 * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
27 * CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
28 * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
29 * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
30 * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
31 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
32 * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
33 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
34 */
35
36 #include <math.h>
37
38 #include "sweep.h"
39 #include "sweep_context.h"
40 #include "advancing_front.h"
41 #include "../common/utils.h"
42 #include "../common/shapes.h"
43
44 void
45 p2t_sweep_init (P2tSweep* THIS)
46 {
47 THIS->nodes_ = g_ptr_array_new ();
48 }
49
50 P2tSweep*
51 p2t_sweep_new ()
52 {
53 P2tSweep* THIS = g_new (P2tSweep, 1);
54 p2t_sweep_init (THIS);
55 return THIS;
56 }
57
58 /**
59 * Destructor - clean up memory
60 */
61 void
62 p2t_sweep_destroy (P2tSweep* THIS)
63 {
64 int i;
65 /* Clean up memory */
66 for (i = 0; i < THIS->nodes_->len; i++)
67 {
68 p2t_node_free (node_index (THIS->nodes_, i));
69 }
70
71 g_ptr_array_free (THIS->nodes_, TRUE);
72 }
73
74 void
75 p2t_sweep_free (P2tSweep* THIS)
76 {
77 p2t_sweep_destroy (THIS);
78 g_free (THIS);
79 }
80
81 /* Triangulate simple polygon with holes */
82
83 void
84 p2t_sweep_triangulate (P2tSweep *THIS, P2tSweepContext *tcx)
85 {
86 p2t_sweepcontext_init_triangulation (tcx);
87 p2t_sweepcontext_create_advancingfront (tcx, THIS->nodes_);
88 /* Sweep points; build mesh */
89 p2t_sweep_sweep_points (THIS, tcx);
90 /* Clean up */
91 p2t_sweep_finalization_polygon (THIS, tcx);
92 }
93
94 void
95 p2t_sweep_sweep_points (P2tSweep *THIS, P2tSweepContext *tcx)
96 {
97 stack_pt++;
98 if (stack_pt > STACKPTMAX)
99 {
100 THROW( MAPTOOL_00001_EXCEPTION );
101 }
102
103 int i, j;
104 for (i = 1; i < p2t_sweepcontext_point_count (tcx); i++)
105 {
106 P2tPoint* point = p2t_sweepcontext_get_point (tcx, i);
107 P2tNode* node = p2t_sweep_point_event (THIS, tcx, point);
108 for (j = 0; j < point->edge_list->len; j++)
109 {
110 p2t_sweep_edge_event_ed_n (THIS, tcx, edge_index (point->edge_list, j), node);
111 }
112 }
113 stack_pt--;
114 }
115
116 void
117 p2t_sweep_finalization_polygon (P2tSweep *THIS, P2tSweepContext *tcx)
118 {
119 /* Get an Internal triangle to start with */
120 P2tTriangle* t = p2t_advancingfront_head (p2t_sweepcontext_front (tcx))->next->triangle;
121 P2tPoint* p = p2t_advancingfront_head (p2t_sweepcontext_front (tcx))->next->point;
122 while (!p2t_triangle_get_constrained_edge_cw (t, p))
123 {
124 t = p2t_triangle_neighbor_ccw (t, p);
125 }
126
127 /* Collect interior triangles constrained by edges */
128 p2t_sweepcontext_mesh_clean (tcx, t);
129 }
130
131 P2tNode*
132 p2t_sweep_point_event (P2tSweep *THIS, P2tSweepContext *tcx, P2tPoint* point)
133 {
134 P2tNode* node = p2t_sweepcontext_locate_node (tcx, point);
135 P2tNode* new_node = p2t_sweep_new_front_triangle (THIS, tcx, point, node);
136
137 /* Only need to check +epsilon since point never have smaller
138 * x value than node due to how we fetch nodes from the front */
139 if (point->x <= node->point->x + EPSILON)
140 {
141 p2t_sweep_fill (THIS, tcx, node);
142 }
143
144 /*tcx.AddNode(new_node); */
145
146 p2t_sweep_fill_advancingfront (THIS, tcx, new_node);
147 return new_node;
148 }
149
150 void
151 p2t_sweep_edge_event_ed_n (P2tSweep *THIS, P2tSweepContext *tcx, P2tEdge* edge, P2tNode* node)
152 {
153 tcx->edge_event.constrained_edge = edge;
154 tcx->edge_event.right = (edge->p->x > edge->q->x);
155
156 if (node == NULL)
157 {
158 //printf("EE005a\n");
159 THROW( MAPTOOL_00001_EXCEPTION );
160 }
161
162 if (node->triangle == NULL)
163 {
164 //printf("EE005b\n");
165 THROW( MAPTOOL_00001_EXCEPTION );
166 }
167
168
169 if (p2t_sweep_is_edge_side_of_triangle (THIS, node->triangle, edge->p, edge->q))
170 {
171 return;
172 }
173
174 /* For now we will do all needed filling
175 * TODO: integrate with flip process might give some better performance
176 * but for now this avoid the issue with cases that needs both flips and fills
177 */
178
179
180 p2t_sweep_fill_edge_event (THIS, tcx, edge, node);
181 p2t_sweep_edge_event_pt_pt_tr_pt (THIS, tcx, edge->p, edge->q, node->triangle, edge->q);
182 }
183
184 void
185 p2t_sweep_edge_event_pt_pt_tr_pt (P2tSweep *THIS, P2tSweepContext *tcx, P2tPoint* ep, P2tPoint* eq, P2tTriangle* triangle, P2tPoint* point)
186 {
187 P2tPoint *p1, *p2;
188 P2tOrientation o1, o2;
189
190 if (triangle == NULL)
191 {
192 //printf("EE001\n");
193 THROW( MAPTOOL_00001_EXCEPTION );
194 }
195
196 if (p2t_sweep_is_edge_side_of_triangle (THIS, triangle, ep, eq))
197 {
198 return;
199 }
200
201 p1 = p2t_triangle_point_ccw (triangle, point);
202 o1 = p2t_orient2d (eq, p1, ep);
203 if (o1 == COLLINEAR)
204 {
205 if (p2t_triangle_contains_pt_pt (triangle, eq, p1))
206 {
207 p2t_triangle_mark_constrained_edge_pt_pt (triangle, eq, p1);
208 /* We are modifying the constraint maybe it would be better to
209 * not change the given constraint and just keep a variable for the new constraint
210 */
211 tcx->edge_event.constrained_edge->q = p1;
212 triangle = p2t_triangle_neighbor_across (triangle, point);
213 p2t_sweep_edge_event_pt_pt_tr_pt (THIS, tcx, ep, p1, triangle, p1);
214 }
215 else
216 {
217 g_error ("EdgeEvent - collinear points not supported");
218 }
219 return;
220 }
221
222 p2 = p2t_triangle_point_cw (triangle, point);
223 o2 = p2t_orient2d (eq, p2, ep);
224 if (o2 == COLLINEAR)
225 {
226 if (p2t_triangle_contains_pt_pt (triangle, eq, p2))
227 {
228 p2t_triangle_mark_constrained_edge_pt_pt (triangle, eq, p2);
229 /* We are modifying the constraint maybe it would be better to
230 * not change the given constraint and just keep a variable for the new constraint
231 */
232 tcx->edge_event.constrained_edge->q = p2;
233 triangle = p2t_triangle_neighbor_across (triangle, point);
234 p2t_sweep_edge_event_pt_pt_tr_pt (THIS, tcx, ep, p2, triangle, p2);
235 }
236 else
237 {
238 g_error ("EdgeEvent - collinear points not supported");
239 }
240 return;
241 }
242
243 if (o1 == o2)
244 {
245 /* Need to decide if we are rotating CW or CCW to get to a triangle
246 * that will cross edge */
247 if (o1 == CW)
248 {
249 triangle = p2t_triangle_neighbor_ccw (triangle, point);
250 }
251 else
252 {
253 triangle = p2t_triangle_neighbor_cw (triangle, point);
254 }
255 p2t_sweep_edge_event_pt_pt_tr_pt (THIS, tcx, ep, eq, triangle, point);
256 }
257 else
258 {
259 /* This triangle crosses constraint so lets flippin start! */
260 p2t_sweep_flip_edge_event (THIS, tcx, ep, eq, triangle, point);
261 }
262 }
263
264 gboolean
265 p2t_sweep_is_edge_side_of_triangle (P2tSweep *THIS, P2tTriangle *triangle, P2tPoint* ep, P2tPoint* eq)
266 {
267 int index = p2t_triangle_edge_index (triangle, ep, eq);
268
269 if (index != -1)
270 {
271 P2tTriangle *t;
272 p2t_triangle_mark_constrained_edge_i (triangle, index);
273 t = p2t_triangle_get_neighbor (triangle, index);
274 if (t)
275 {
276 p2t_triangle_mark_constrained_edge_pt_pt (t, ep, eq);
277 }
278 return TRUE;
279 }
280 return FALSE;
281 }
282
283 P2tNode*
284 p2t_sweep_new_front_triangle (P2tSweep *THIS, P2tSweepContext *tcx, P2tPoint* point, P2tNode *node)
285 {
286 P2tTriangle* triangle = p2t_triangle_new (point, node->point, node->next->point);
287 P2tNode *new_node;
288
289 p2t_triangle_mark_neighbor_tr (triangle, node->triangle);
290 p2t_sweepcontext_add_to_map (tcx, triangle);
291
292 new_node = p2t_node_new_pt (point);
293 g_ptr_array_add (THIS->nodes_, new_node);
294
295 new_node->next = node->next;
296 new_node->prev = node;
297 node->next->prev = new_node;
298 node->next = new_node;
299
300 if (!p2t_sweep_legalize (THIS, tcx, triangle))
301 {
302 p2t_sweepcontext_map_triangle_to_nodes (tcx, triangle);
303 }
304
305 return new_node;
306 }
307
308 void
309 p2t_sweep_fill (P2tSweep *THIS, P2tSweepContext *tcx, P2tNode* node)
310 {
311 P2tTriangle* triangle = p2t_triangle_new (node->prev->point, node->point, node->next->point);
312
313 /* TODO: should copy the constrained_edge value from neighbor triangles
314 * for now constrained_edge values are copied during the legalize */
315 p2t_triangle_mark_neighbor_tr (triangle, node->prev->triangle);
316 p2t_triangle_mark_neighbor_tr (triangle, node->triangle);
317
318 p2t_sweepcontext_add_to_map (tcx, triangle);
319
320 /* Update the advancing front */
321 node->prev->next = node->next;
322 node->next->prev = node->prev;
323
324 /* If it was legalized the triangle has already been mapped */
325 if (!p2t_sweep_legalize (THIS, tcx, triangle))
326 {
327 p2t_sweepcontext_map_triangle_to_nodes (tcx, triangle);
328 }
329
330 }
331
332 void
333 p2t_sweep_fill_advancingfront (P2tSweep *THIS, P2tSweepContext *tcx, P2tNode* n)
334 {
335
336 /* Fill right holes */
337 P2tNode* node = n->next;
338
339 while (node->next)
340 {
341 /* if HoleAngle exceeds 90 degrees then break. */
342 if (p2t_sweep_large_hole_dont_fill (THIS, node)) break;
343 p2t_sweep_fill (THIS, tcx, node);
344 node = node->next;
345 }
346
347 /* Fill left holes */
348 node = n->prev;
349
350 while (node->prev)
351 {
352 /* if HoleAngle exceeds 90 degrees then break. */
353 if (p2t_sweep_large_hole_dont_fill (THIS, node)) break;
354 p2t_sweep_fill (THIS, tcx, node);
355 node = node->prev;
356 }
357
358 /* Fill right basins */
359 if (n->next && n->next->next)
360 {
361 double angle = p2t_sweep_basin_angle (THIS, n);
362 if (angle < PI_3div4)
363 {
364 p2t_sweep_fill_basin (THIS, tcx, n);
365 }
366 }
367 }
368
369 /* True if HoleAngle exceeds 90 degrees. */
370 gboolean
371 p2t_sweep_large_hole_dont_fill (P2tSweep *THIS, P2tNode* node)
372 {
373 P2tNode* nextNode = node->next;
374 P2tNode* prevNode = node->prev;
375 P2tNode *next2Node, *prev2Node;
376 if (! p2t_sweep_angle_exceeds_90_degrees (THIS, node->point, nextNode->point, prevNode->point))
377 return FALSE;
378
379 /* Check additional points on front. */
380 next2Node = nextNode->next;
381 /* "..Plus.." because only want angles on same side as point being added. */
382 if ((next2Node != NULL) && !p2t_sweep_angle_exceeds_plus_90_degrees_or_is_negative (THIS, node->point, next2Node->point, prevNode->point))
383 return FALSE;
384
385 prev2Node = prevNode->prev;
386 /* "..Plus.." because only want angles on same side as point being added. */
387 if ((prev2Node != NULL) && !p2t_sweep_angle_exceeds_plus_90_degrees_or_is_negative (THIS, node->point, nextNode->point, prev2Node->point))
388 return FALSE;
389
390 return TRUE;
391 }
392
393 gboolean
394 p2t_sweep_angle_exceeds_90_degrees(P2tSweep* THIS, P2tPoint* origin, P2tPoint* pa, P2tPoint* pb)
395 {
396 gdouble angle = p2t_sweep_angle (THIS, origin, pa, pb);
397 gboolean exceeds90Degrees = ((angle > G_PI_2) || (angle < -G_PI_2));
398 return exceeds90Degrees;
399 }
400
401 gboolean
402 p2t_sweep_angle_exceeds_plus_90_degrees_or_is_negative (P2tSweep* THIS, P2tPoint* origin, P2tPoint* pa, P2tPoint* pb)
403 {
404 gdouble angle = p2t_sweep_angle (THIS, origin, pa, pb);
405 gboolean exceedsPlus90DegreesOrIsNegative = (angle > G_PI_2) || (angle < 0);
406 return exceedsPlus90DegreesOrIsNegative;
407 }
408
409 gdouble
410 p2t_sweep_angle (P2tSweep* THIS, P2tPoint* origin, P2tPoint* pa, P2tPoint* pb) {
411 /* Complex plane
412 * ab = cosA +i*sinA
413 * ab = (ax + ay*i)(bx + by*i) = (ax*bx + ay*by) + i(ax*by-ay*bx)
414 * atan2(y,x) computes the principal value of the argument function
415 * applied to the complex number x+iy
416 * Where x = ax*bx + ay*by
417 * y = ax*by - ay*bx
418 */
419 double px = origin->x;
420 double py = origin->y;
421 double ax = pa->x - px;
422 double ay = pa->y - py;
423 double bx = pb->x - px;
424 double by = pb->y - py;
425 double x = ax * by - ay * bx;
426 double y = ax * bx + ay * by;
427 double angle = atan2(x, y);
428 return angle;
429 }
430
431 double
432 p2t_sweep_basin_angle (P2tSweep *THIS, P2tNode* node)
433 {
434 double ax = node->point->x - node->next->next->point->x;
435 double ay = node->point->y - node->next->next->point->y;
436 return atan2 (ay, ax);
437 }
438
439 double
440 p2t_sweep_hole_angle (P2tSweep *THIS, P2tNode* node)
441 {
442 /* Complex plane
443 * ab = cosA +i*sinA
444 * ab = (ax + ay*i)(bx + by*i) = (ax*bx + ay*by) + i(ax*by-ay*bx)
445 * atan2(y,x) computes the principal value of the argument function
446 * applied to the complex number x+iy
447 * Where x = ax*bx + ay*by
448 * y = ax*by - ay*bx
449 */
450 double ax = node->next->point->x - node->point->x;
451 double ay = node->next->point->y - node->point->y;
452 double bx = node->prev->point->x - node->point->x;
453 double by = node->prev->point->y - node->point->y;
454 return atan2 (ax * by - ay * bx, ax * bx + ay * by);
455 }
456
457 gboolean
458 p2t_sweep_legalize (P2tSweep *THIS, P2tSweepContext *tcx, P2tTriangle *t)
459 {
460 int i;
461 /* To legalize a triangle we start by finding if any of the three edges
462 * violate the Delaunay condition */
463 for (i = 0; i < 3; i++)
464 {
465 P2tTriangle *ot;
466
467 if (t->delaunay_edge[i])
468 continue;
469
470 ot = p2t_triangle_get_neighbor (t, i);
471
472 if (ot)
473 {
474 P2tPoint* p = p2t_triangle_get_point (t, i);
475 P2tPoint* op = p2t_triangle_opposite_point (ot, t, p);
476 int oi = p2t_triangle_index (ot, op);
477 gboolean inside;
478
479 /* If this is a Constrained Edge or a Delaunay Edge(only during recursive legalization)
480 * then we should not try to legalize */
481 if (ot->constrained_edge[oi] || ot->delaunay_edge[oi])
482 {
483 t->constrained_edge[i] = ot->constrained_edge[oi];
484 continue;
485 }
486
487 inside = p2t_sweep_incircle (THIS, p, p2t_triangle_point_ccw (t, p), p2t_triangle_point_cw (t, p), op);
488
489 if (inside)
490 {
491 gboolean not_legalized;
492 /* Lets mark this shared edge as Delaunay */
493 t->delaunay_edge[i] = TRUE;
494 ot->delaunay_edge[oi] = TRUE;
495
496 /* Lets rotate shared edge one vertex CW to legalize it */
497 p2t_sweep_rotate_triangle_pair (THIS, t, p, ot, op);
498
499 /* We now got one valid Delaunay Edge shared by two triangles
500 * This gives us 4 new edges to check for Delaunay */
501
502 /* Make sure that triangle to node mapping is done only one time for a specific triangle */
503 not_legalized = !p2t_sweep_legalize (THIS, tcx, t);
504 if (not_legalized)
505 {
506 p2t_sweepcontext_map_triangle_to_nodes (tcx, t);
507 }
508
509 not_legalized = !p2t_sweep_legalize (THIS, tcx, ot);
510 if (not_legalized)
511 p2t_sweepcontext_map_triangle_to_nodes (tcx, ot);
512
513 /* Reset the Delaunay edges, since they only are valid Delaunay edges
514 * until we add a new triangle or point.
515 * XXX: need to think about this. Can these edges be tried after we
516 * return to previous recursive level? */
517 t->delaunay_edge[i] = FALSE;
518 ot->delaunay_edge[oi] = FALSE;
519
520 /* If triangle have been legalized no need to check the other edges since
521 * the recursive legalization will handles those so we can end here.*/
522 return TRUE;
523 }
524 }
525 }
526 return FALSE;
527 }
528
529 gboolean
530 p2t_sweep_incircle (P2tSweep *THIS, P2tPoint* pa, P2tPoint* pb, P2tPoint* pc, P2tPoint* pd)
531 {
532 double adx = pa->x - pd->x;
533 double ady = pa->y - pd->y;
534 double bdx = pb->x - pd->x;
535 double bdy = pb->y - pd->y;
536
537 double adxbdy = adx * bdy;
538 double bdxady = bdx * ady;
539 double oabd = adxbdy - bdxady;
540
541 double cdx, cdy;
542 double cdxady, adxcdy, ocad;
543
544 double bdxcdy, cdxbdy;
545 double alift, blift, clift;
546 double det;
547
548 if (oabd <= 0)
549 return FALSE;
550
551 cdx = pc->x - pd->x;
552 cdy = pc->y - pd->y;
553
554 cdxady = cdx * ady;
555 adxcdy = adx * cdy;
556 ocad = cdxady - adxcdy;
557
558 if (ocad <= 0)
559 return FALSE;
560
561 bdxcdy = bdx * cdy;
562 cdxbdy = cdx * bdy;
563
564 alift = adx * adx + ady * ady;
565 blift = bdx * bdx + bdy * bdy;
566 clift = cdx * cdx + cdy * cdy;
567
568 det = alift * (bdxcdy - cdxbdy) + blift * ocad + clift * oabd;
569
570 return det > 0;
571 }
572
573 void
574 p2t_sweep_rotate_triangle_pair (P2tSweep *THIS, P2tTriangle *t, P2tPoint* p, P2tTriangle *ot, P2tPoint* op)
575 {
576 P2tTriangle *n1, *n2, *n3, *n4;
577 gboolean ce1, ce2, ce3, ce4;
578 gboolean de1, de2, de3, de4;
579
580 n1 = p2t_triangle_neighbor_ccw (t, p);
581 n2 = p2t_triangle_neighbor_cw (t, p);
582 n3 = p2t_triangle_neighbor_ccw (ot, op);
583 n4 = p2t_triangle_neighbor_cw (ot, op);
584
585 ce1 = p2t_triangle_get_constrained_edge_ccw (t, p);
586 ce2 = p2t_triangle_get_constrained_edge_cw (t, p);
587 ce3 = p2t_triangle_get_constrained_edge_ccw (ot, op);
588 ce4 = p2t_triangle_get_constrained_edge_cw (ot, op);
589
590 de1 = p2t_triangle_get_delunay_edge_ccw (t, p);
591 de2 = p2t_triangle_get_delunay_edge_cw (t, p);
592 de3 = p2t_triangle_get_delunay_edge_ccw (ot, op);
593 de4 = p2t_triangle_get_delunay_edge_cw (ot, op);
594
595 p2t_triangle_legalize_pt_pt (t, p, op);
596 p2t_triangle_legalize_pt_pt (ot, op, p);
597
598 /* Remap delaunay_edge */
599 p2t_triangle_set_delunay_edge_ccw (ot, p, de1);
600 p2t_triangle_set_delunay_edge_cw (t, p, de2);
601 p2t_triangle_set_delunay_edge_ccw (t, op, de3);
602 p2t_triangle_set_delunay_edge_cw (ot, op, de4);
603
604 /* Remap constrained_edge */
605 p2t_triangle_set_constrained_edge_ccw (ot, p, ce1);
606 p2t_triangle_set_constrained_edge_cw (t, p, ce2);
607 p2t_triangle_set_constrained_edge_ccw (t, op, ce3);
608 p2t_triangle_set_constrained_edge_cw (ot, op, ce4);
609
610 /* Remap neighbors
611 * XXX: might optimize the markNeighbor by keeping track of
612 * what side should be assigned to what neighbor after the
613 * rotation. Now mark neighbor does lots of testing to find
614 * the right side. */
615 p2t_triangle_clear_neighbors (t);
616 p2t_triangle_clear_neighbors (ot);
617 if (n1) p2t_triangle_mark_neighbor_tr (ot, n1);
618 if (n2) p2t_triangle_mark_neighbor_tr (t, n2);
619 if (n3) p2t_triangle_mark_neighbor_tr (t, n3);
620 if (n4) p2t_triangle_mark_neighbor_tr (ot, n4);
621 p2t_triangle_mark_neighbor_tr (t, ot);
622 }
623
624 void
625 p2t_sweep_fill_basin (P2tSweep *THIS, P2tSweepContext *tcx, P2tNode* node)
626 {
627 if (p2t_orient2d (node->point, node->next->point, node->next->next->point) == CCW)
628 {
629 tcx->basin.left_node = node->next->next;
630 }
631 else
632 {
633 tcx->basin.left_node = node->next;
634 }
635
636 /* Find the bottom and right node */
637 tcx->basin.bottom_node = tcx->basin.left_node;
638 while (tcx->basin.bottom_node->next
639 && tcx->basin.bottom_node->point->y >= tcx->basin.bottom_node->next->point->y)
640 {
641 tcx->basin.bottom_node = tcx->basin.bottom_node->next;
642 }
643 if (tcx->basin.bottom_node == tcx->basin.left_node)
644 {
645 /* No valid basin */
646 return;
647 }
648
649 tcx->basin.right_node = tcx->basin.bottom_node;
650 while (tcx->basin.right_node->next
651 && tcx->basin.right_node->point->y < tcx->basin.right_node->next->point->y)
652 {
653 tcx->basin.right_node = tcx->basin.right_node->next;
654 }
655 if (tcx->basin.right_node == tcx->basin.bottom_node)
656 {
657 /* No valid basins */
658 return;
659 }
660
661 tcx->basin.width = tcx->basin.right_node->point->x - tcx->basin.left_node->point->x;
662 tcx->basin.left_highest = tcx->basin.left_node->point->y > tcx->basin.right_node->point->y;
663
664 p2t_sweep_fill_basin_req (THIS, tcx, tcx->basin.bottom_node);
665 }
666
667 void
668 p2t_sweep_fill_basin_req (P2tSweep *THIS, P2tSweepContext *tcx, P2tNode* node)
669 {
670 /* if shallow stop filling */
671 if (p2t_sweep_is_shallow (THIS, tcx, node))
672 {
673 return;
674 }
675
676 p2t_sweep_fill (THIS, tcx, node);
677
678 if (node->prev == tcx->basin.left_node && node->next == tcx->basin.right_node)
679 {
680 return;
681 }
682 else if (node->prev == tcx->basin.left_node)
683 {
684 P2tOrientation o = p2t_orient2d (node->point, node->next->point, node->next->next->point);
685 if (o == CW)
686 {
687 return;
688 }
689 node = node->next;
690 }
691 else if (node->next == tcx->basin.right_node)
692 {
693 P2tOrientation o = p2t_orient2d (node->point, node->prev->point, node->prev->prev->point);
694 if (o == CCW)
695 {
696 return;
697 }
698 node = node->prev;
699 }
700 else
701 {
702 /* Continue with the neighbor node with lowest Y value */
703 if (node->prev->point->y < node->next->point->y)
704 {
705 node = node->prev;
706 }
707 else
708 {
709 node = node->next;
710 }
711 }
712
713 p2t_sweep_fill_basin_req (THIS, tcx, node);
714 }
715
716 gboolean
717 p2t_sweep_is_shallow (P2tSweep *THIS, P2tSweepContext *tcx, P2tNode* node)
718 {
719
720
721 double height;
722
723 if (tcx->basin.left_highest)
724 {
725 height = tcx->basin.left_node->point->y - node->point->y;
726 }
727 else
728 {
729 height = tcx->basin.right_node->point->y - node->point->y;
730 }
731
732 /* if shallow stop filling */
733 if (tcx->basin.width > height)
734 {
735 return TRUE;
736 }
737 return FALSE;
738 }
739
740 void
741 p2t_sweep_fill_edge_event (P2tSweep *THIS, P2tSweepContext *tcx, P2tEdge* edge, P2tNode* node)
742 {
743
744
745 if (tcx->edge_event.right)
746 {
747 p2t_sweep_fill_right_above_edge_event (THIS, tcx, edge, node);
748 }
749 else
750 {
751 p2t_sweep_fill_left_above_edge_event (THIS, tcx, edge, node);
752 }
753 }
754
755 void
756 p2t_sweep_fill_right_above_edge_event (P2tSweep *THIS, P2tSweepContext *tcx, P2tEdge* edge, P2tNode* node)
757 {
758 while (node->next->point->x < edge->p->x)
759 {
760 /* Check if next node is below the edge */
761 if (p2t_orient2d (edge->q, node->next->point, edge->p) == CCW)
762 {
763 p2t_sweep_fill_right_below_edge_event (THIS, tcx, edge, node);
764 }
765 else
766 {
767 node = node->next;
768 }
769 }
770 }
771
772 void
773 p2t_sweep_fill_right_below_edge_event (P2tSweep *THIS, P2tSweepContext *tcx, P2tEdge* edge, P2tNode* node)
774 {
775 if (node->point->x < edge->p->x)
776 {
777 if (p2t_orient2d (node->point, node->next->point, node->next->next->point) == CCW)
778 {
779 /* Concave */
780 p2t_sweep_fill_right_concave_edge_event (THIS, tcx, edge, node);
781 }
782 else
783 {
784 /* Convex */
785 p2t_sweep_fill_right_convex_edge_event (THIS, tcx, edge, node);
786 /* Retry this one */
787 p2t_sweep_fill_right_below_edge_event (THIS, tcx, edge, node);
788 }
789 }
790 }
791
792 void
793 p2t_sweep_fill_right_concave_edge_event (P2tSweep *THIS, P2tSweepContext *tcx, P2tEdge* edge, P2tNode* node)
794 {
795 p2t_sweep_fill (THIS, tcx, node->next);
796 if (node->next->point != edge->p)
797 {
798 /* Next above or below edge? */
799 if (p2t_orient2d (edge->q, node->next->point, edge->p) == CCW)
800 {
801 /* Below */
802 if (p2t_orient2d (node->point, node->next->point, node->next->next->point) == CCW)
803 {
804 /* Next is concave */
805 p2t_sweep_fill_right_concave_edge_event (THIS, tcx, edge, node);
806 }
807 else
808 {
809 /* Next is convex */
810 }
811 }
812 }
813
814 }
815
816 void
817 p2t_sweep_fill_right_convex_edge_event (P2tSweep *THIS, P2tSweepContext *tcx, P2tEdge* edge, P2tNode* node)
818 {
819 /* Next concave or convex? */
820 if (p2t_orient2d (node->next->point, node->next->next->point, node->next->next->next->point) == CCW)
821 {
822 /* Concave */
823 p2t_sweep_fill_right_concave_edge_event (THIS, tcx, edge, node->next);
824 }
825 else
826 {
827 /* Convex
828 * Next above or below edge? */
829 if (p2t_orient2d (edge->q, node->next->next->point, edge->p) == CCW)
830 {
831 /* Below */
832 p2t_sweep_fill_right_convex_edge_event (THIS, tcx, edge, node->next);
833 }
834 else
835 {
836 /* Above */
837 }
838 }
839 }
840
841 void
842 p2t_sweep_fill_left_above_edge_event (P2tSweep *THIS, P2tSweepContext *tcx, P2tEdge* edge, P2tNode* node)
843 {
844 while (node->prev->point->x > edge->p->x)
845 {
846 /* Check if next node is below the edge */
847 if (p2t_orient2d (edge->q, node->prev->point, edge->p) == CW)
848 {
849 p2t_sweep_fill_left_below_edge_event (THIS, tcx, edge, node);
850 }
851 else
852 {
853 node = node->prev;
854 }
855 }
856 }
857
858 void
859 p2t_sweep_fill_left_below_edge_event (P2tSweep *THIS, P2tSweepContext *tcx, P2tEdge* edge, P2tNode* node)
860 {
861 if (node->point->x > edge->p->x)
862 {
863 if (p2t_orient2d (node->point, node->prev->point, node->prev->prev->point) == CW)
864 {
865 /* Concave */
866 p2t_sweep_fill_left_concave_edge_event (THIS, tcx, edge, node);
867 }
868 else
869 {
870 /* Convex */
871 p2t_sweep_fill_left_convex_edge_event (THIS, tcx, edge, node);
872 /* Retry this one */
873 p2t_sweep_fill_left_below_edge_event (THIS, tcx, edge, node);
874 }
875 }
876 }
877
878 void
879 p2t_sweep_fill_left_convex_edge_event (P2tSweep *THIS, P2tSweepContext *tcx, P2tEdge* edge, P2tNode* node)
880 {
881
882
883 /* Next concave or convex? */
884 if (p2t_orient2d (node->prev->point, node->prev->prev->point, node->prev->prev->prev->point) == CW)
885 {
886 /* Concave */
887 p2t_sweep_fill_left_concave_edge_event (THIS, tcx, edge, node->prev);
888 }
889 else
890 {
891 /* Convex
892 * Next above or below edge? */
893 if (p2t_orient2d (edge->q, node->prev->prev->point, edge->p) == CW)
894 {
895 /* Below */
896 p2t_sweep_fill_left_convex_edge_event (THIS, tcx, edge, node->prev);
897 }
898 else
899 {
900 /* Above */
901 }
902 }
903 }
904
905 void
906 p2t_sweep_fill_left_concave_edge_event (P2tSweep *THIS, P2tSweepContext *tcx, P2tEdge* edge, P2tNode* node)
907 {
908
909
910 p2t_sweep_fill (THIS, tcx, node->prev);
911 if (node->prev->point != edge->p)
912 {
913 /* Next above or below edge? */
914 if (p2t_orient2d (edge->q, node->prev->point, edge->p) == CW)
915 {
916 /* Below */
917 if (p2t_orient2d (node->point, node->prev->point, node->prev->prev->point) == CW)
918 {
919 /* Next is concave */
920 p2t_sweep_fill_left_concave_edge_event (THIS, tcx, edge, node);
921 }
922 else
923 {
924 /* Next is convex */
925 }
926 }
927 }
928
929 }
930
931 void
932 p2t_sweep_flip_edge_event (P2tSweep *THIS, P2tSweepContext *tcx, P2tPoint* ep, P2tPoint* eq, P2tTriangle* t, P2tPoint* p)
933 {
934 stack_pt++;
935 if (stack_pt > STACKPTMAX)
936 {
937 THROW( MAPTOOL_00001_EXCEPTION );
938 }
939
940
941 P2tTriangle* ot = p2t_triangle_neighbor_across (t, p);
942 P2tPoint* op = p2t_triangle_opposite_point (ot, t, p);
943
944 if (ot == NULL)
945 {
946 /* If we want to integrate the fillEdgeEvent do it here
947 * With current implementation we should never get here
948 *throw new RuntimeException( "[BUG:FIXME] FLIP failed due to missing triangle");
949 */
950 fprintf(stderr, "EE002\n");
951 THROW( MAPTOOL_00001_EXCEPTION );
952
953 assert (0);
954 }
955
956
957
958 if (p2t_utils_in_scan_area (p, p2t_triangle_point_ccw (t, p), p2t_triangle_point_cw (t, p), op))
959 {
960
961 /* Lets rotate shared edge one vertex CW */
962 p2t_sweep_rotate_triangle_pair (THIS, t, p, ot, op);
963 p2t_sweepcontext_map_triangle_to_nodes (tcx, t);
964 p2t_sweepcontext_map_triangle_to_nodes (tcx, ot);
965
966
967 if (p == eq && op == ep)
968 {
969 if (p2t_point_equals (eq, tcx->edge_event.constrained_edge->q) && p2t_point_equals (ep, tcx->edge_event.constrained_edge->p))
970 {
971
972 p2t_triangle_mark_constrained_edge_pt_pt (t, ep, eq);
973 p2t_triangle_mark_constrained_edge_pt_pt (ot, ep, eq);
974 p2t_sweep_legalize (THIS, tcx, t);
975 p2t_sweep_legalize (THIS, tcx, ot);
976
977 }
978 else
979 {
980 /* XXX: I think one of the triangles should be legalized here? */
981
982 }
983 }
984 else
985 {
986
987 P2tOrientation o = p2t_orient2d (eq, op, ep);
988 t = p2t_sweep_next_flip_triangle (THIS, tcx, (int) o, t, ot, p, op);
989 p2t_sweep_flip_edge_event (THIS, tcx, ep, eq, t, p);
990
991 }
992 }
993 else
994 {
995
996 P2tPoint* newP = p2t_sweep_next_flip_point (THIS, ep, eq, ot, op);
997
998 p2t_sweep_flip_scan_edge_event (THIS, tcx, ep, eq, t, ot, newP);
999
1000 p2t_sweep_edge_event_pt_pt_tr_pt (THIS, tcx, ep, eq, t, p);
1001
1002 }
1003
1004 stack_pt--;
1005
1006 }
1007
1008 P2tTriangle*
1009 p2t_sweep_next_flip_triangle (P2tSweep *THIS, P2tSweepContext *tcx, int o, P2tTriangle *t, P2tTriangle *ot, P2tPoint* p, P2tPoint* op)
1010 {
1011
1012
1013 int edge_index;
1014
1015 if (o == CCW)
1016 {
1017 /* ot is not crossing edge after flip */
1018 int edge_index = p2t_triangle_edge_index (ot, p, op);
1019 ot->delaunay_edge[edge_index] = TRUE;
1020 p2t_sweep_legalize (THIS, tcx, ot);
1021 p2t_triangle_clear_delunay_edges (ot);
1022 return t;
1023 }
1024
1025 /* t is not crossing edge after flip */
1026 edge_index = p2t_triangle_edge_index (t, p, op);
1027
1028 t->delaunay_edge[edge_index] = TRUE;
1029 p2t_sweep_legalize (THIS, tcx, t);
1030 p2t_triangle_clear_delunay_edges (t);
1031 return ot;
1032 }
1033
1034 P2tPoint*
1035 p2t_sweep_next_flip_point (P2tSweep *THIS, P2tPoint* ep, P2tPoint* eq, P2tTriangle *ot, P2tPoint* op)
1036 {
1037
1038 stack_pt++;
1039 if (stack_pt > STACKPTMAX)
1040 {
1041 THROW( MAPTOOL_00001_EXCEPTION );
1042 }
1043
1044
1045 P2tOrientation o2d = p2t_orient2d (eq, op, ep);
1046 if (o2d == CW)
1047 {
1048 /* Right */
1049 return p2t_triangle_point_ccw (ot, op);
1050 }
1051 else if (o2d == CCW)
1052 {
1053 /* Left */
1054 return p2t_triangle_point_cw (ot, op);
1055 }
1056 else
1057 {
1058 /*throw new RuntimeException("[Unsupported] Opposing point on constrained edge");*/
1059 //printf("EE003\n");
1060 THROW( MAPTOOL_00001_EXCEPTION );
1061
1062 assert (0);
1063 }
1064
1065 stack_pt--;
1066 }
1067
1068 void
1069 p2t_sweep_flip_scan_edge_event (P2tSweep *THIS, P2tSweepContext *tcx, P2tPoint* ep, P2tPoint* eq, P2tTriangle *flip_triangle,
1070 P2tTriangle *t, P2tPoint* p)
1071 {
1072
1073 stack_pt++;
1074 if (stack_pt > STACKPTMAX)
1075 {
1076 THROW( MAPTOOL_00001_EXCEPTION );
1077 }
1078
1079
1080 P2tTriangle* ot = p2t_triangle_neighbor_across (t, p);
1081
1082 P2tPoint* op = p2t_triangle_opposite_point (ot, t, p);
1083
1084
1085 if (p2t_triangle_neighbor_across (t, p) == NULL)
1086 {
1087 /* If we want to integrate the fillEdgeEvent do it here
1088 * With current implementation we should never get here
1089 *throw new RuntimeException( "[BUG:FIXME] FLIP failed due to missing triangle");
1090 */
1091 //printf("EE004\n");
1092 THROW( MAPTOOL_00001_EXCEPTION );
1093
1094 assert (0);
1095 }
1096
1097
1098 if (p2t_utils_in_scan_area (eq, p2t_triangle_point_ccw (flip_triangle, eq), p2t_triangle_point_cw (flip_triangle, eq), op))
1099 {
1100 /* flip with new edge op->eq */
1101
1102 p2t_sweep_flip_edge_event (THIS, tcx, eq, op, ot, op);
1103
1104 /* TODO: Actually I just figured out that it should be possible to
1105 * improve this by getting the next ot and op before the the above
1106 * flip and continue the flipScanEdgeEvent here
1107 * set new ot and op here and loop back to inScanArea test
1108 * also need to set a new flip_triangle first
1109 * Turns out at first glance that this is somewhat complicated
1110 * so it will have to wait. */
1111 }
1112 else
1113 {
1114
1115 P2tPoint* newP = p2t_sweep_next_flip_point (THIS, ep, eq, ot, op);
1116
1117 p2t_sweep_flip_scan_edge_event (THIS, tcx, ep, eq, flip_triangle, ot, newP);
1118
1119 }
1120
1121 stack_pt--;
1122
1123
1124 }
1125

   
Visit the ZANavi Wiki