/[zanavi_public1]/navit/navit/maptool/poly2tri-c/001/seidel-1.0/construct.c
ZANavi

Contents of /navit/navit/maptool/poly2tri-c/001/seidel-1.0/construct.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 31 - (show annotations) (download)
Mon Feb 4 17:41:59 2013 UTC (11 years, 1 month ago) by zoff99
File MIME type: text/plain
File size: 25069 byte(s)
new map version, lots of fixes and experimental new features
1 #include <triangulate.h>
2 #include <math.h>
3
4
5 node_t qs[QSIZE]; /* Query structure */
6 trap_t tr[TRSIZE]; /* Trapezoid structure */
7 segment_t seg[SEGSIZE]; /* Segment table */
8
9 static int q_idx;
10 static int tr_idx;
11
12 /* Return a new node to be added into the query tree */
13 static int newnode()
14 {
15 if (q_idx < QSIZE)
16 return q_idx++;
17 else
18 {
19 fprintf(stderr, "newnode: Query-table overflow\n");
20 return -1;
21 }
22 }
23
24 /* Return a free trapezoid */
25 static int newtrap()
26 {
27 if (tr_idx < TRSIZE)
28 {
29 tr[tr_idx].lseg = -1;
30 tr[tr_idx].rseg = -1;
31 tr[tr_idx].state = ST_VALID;
32 return tr_idx++;
33 }
34 else
35 {
36 fprintf(stderr, "newtrap: Trapezoid-table overflow\n");
37 return -1;
38 }
39 }
40
41
42 /* Return the maximum of the two points into the yval structure */
43 static int _max(yval, v0, v1)
44 point_t *yval;
45 point_t *v0;
46 point_t *v1;
47 {
48 if (v0->y > v1->y + C_EPS)
49 *yval = *v0;
50 else if (FP_EQUAL(v0->y, v1->y))
51 {
52 if (v0->x > v1->x + C_EPS)
53 *yval = *v0;
54 else
55 *yval = *v1;
56 }
57 else
58 *yval = *v1;
59
60 return 0;
61 }
62
63
64 /* Return the minimum of the two points into the yval structure */
65 static int _min(yval, v0, v1)
66 point_t *yval;
67 point_t *v0;
68 point_t *v1;
69 {
70 if (v0->y < v1->y - C_EPS)
71 *yval = *v0;
72 else if (FP_EQUAL(v0->y, v1->y))
73 {
74 if (v0->x < v1->x)
75 *yval = *v0;
76 else
77 *yval = *v1;
78 }
79 else
80 *yval = *v1;
81
82 return 0;
83 }
84
85
86 int _greater_than(v0, v1)
87 point_t *v0;
88 point_t *v1;
89 {
90 if (v0->y > v1->y + C_EPS)
91 return TRUE;
92 else if (v0->y < v1->y - C_EPS)
93 return FALSE;
94 else
95 return (v0->x > v1->x);
96 }
97
98
99 int _equal_to(v0, v1)
100 point_t *v0;
101 point_t *v1;
102 {
103 return (FP_EQUAL(v0->y, v1->y) && FP_EQUAL(v0->x, v1->x));
104 }
105
106 int _greater_than_equal_to(v0, v1)
107 point_t *v0;
108 point_t *v1;
109 {
110 if (v0->y > v1->y + C_EPS)
111 return TRUE;
112 else if (v0->y < v1->y - C_EPS)
113 return FALSE;
114 else
115 return (v0->x >= v1->x);
116 }
117
118 int _less_than(v0, v1)
119 point_t *v0;
120 point_t *v1;
121 {
122 if (v0->y < v1->y - C_EPS)
123 return TRUE;
124 else if (v0->y > v1->y + C_EPS)
125 return FALSE;
126 else
127 return (v0->x < v1->x);
128 }
129
130
131 /* Initilialise the query structure (Q) and the trapezoid table (T)
132 * when the first segment is added to start the trapezoidation. The
133 * query-tree starts out with 4 trapezoids, one S-node and 2 Y-nodes
134 *
135 * 4
136 * -----------------------------------
137 * \
138 * 1 \ 2
139 * \
140 * -----------------------------------
141 * 3
142 */
143
144 static int init_query_structure(segnum)
145 int segnum;
146 {
147 int i1, i2, i3, i4, i5, i6, i7, root;
148 int t1, t2, t3, t4;
149 segment_t *s = &seg[segnum];
150
151 q_idx = tr_idx = 1;
152 memset((void *)tr, 0, sizeof(tr));
153 memset((void *)qs, 0, sizeof(qs));
154
155 i1 = newnode();
156 qs[i1].nodetype = T_Y;
157 _max(&qs[i1].yval, &s->v0, &s->v1); /* root */
158 root = i1;
159
160 qs[i1].right = i2 = newnode();
161 qs[i2].nodetype = T_SINK;
162 qs[i2].parent = i1;
163
164 qs[i1].left = i3 = newnode();
165 qs[i3].nodetype = T_Y;
166 _min(&qs[i3].yval, &s->v0, &s->v1); /* root */
167 qs[i3].parent = i1;
168
169 qs[i3].left = i4 = newnode();
170 qs[i4].nodetype = T_SINK;
171 qs[i4].parent = i3;
172
173 qs[i3].right = i5 = newnode();
174 qs[i5].nodetype = T_X;
175 qs[i5].segnum = segnum;
176 qs[i5].parent = i3;
177
178 qs[i5].left = i6 = newnode();
179 qs[i6].nodetype = T_SINK;
180 qs[i6].parent = i5;
181
182 qs[i5].right = i7 = newnode();
183 qs[i7].nodetype = T_SINK;
184 qs[i7].parent = i5;
185
186 t1 = newtrap(); /* middle left */
187 t2 = newtrap(); /* middle right */
188 t3 = newtrap(); /* bottom-most */
189 t4 = newtrap(); /* topmost */
190
191 tr[t1].hi = tr[t2].hi = tr[t4].lo = qs[i1].yval;
192 tr[t1].lo = tr[t2].lo = tr[t3].hi = qs[i3].yval;
193 tr[t4].hi.y = (double) (INFINITY);
194 tr[t4].hi.x = (double) (INFINITY);
195 tr[t3].lo.y = (double) -1* (INFINITY);
196 tr[t3].lo.x = (double) -1* (INFINITY);
197 tr[t1].rseg = tr[t2].lseg = segnum;
198 tr[t1].u0 = tr[t2].u0 = t4;
199 tr[t1].d0 = tr[t2].d0 = t3;
200 tr[t4].d0 = tr[t3].u0 = t1;
201 tr[t4].d1 = tr[t3].u1 = t2;
202
203 tr[t1].sink = i6;
204 tr[t2].sink = i7;
205 tr[t3].sink = i4;
206 tr[t4].sink = i2;
207
208 tr[t1].state = tr[t2].state = ST_VALID;
209 tr[t3].state = tr[t4].state = ST_VALID;
210
211 qs[i2].trnum = t4;
212 qs[i4].trnum = t3;
213 qs[i6].trnum = t1;
214 qs[i7].trnum = t2;
215
216 s->is_inserted = TRUE;
217 return root;
218 }
219
220
221 /* Retun TRUE if the vertex v is to the left of line segment no.
222 * segnum. Takes care of the degenerate cases when both the vertices
223 * have the same y--cood, etc.
224 */
225
226 static int is_left_of(segnum, v)
227 int segnum;
228 point_t *v;
229 {
230 segment_t *s = &seg[segnum];
231 double area;
232
233 if (_greater_than(&s->v1, &s->v0)) /* seg. going upwards */
234 {
235 if (FP_EQUAL(s->v1.y, v->y))
236 {
237 if (v->x < s->v1.x)
238 area = 1.0;
239 else
240 area = -1.0;
241 }
242 else if (FP_EQUAL(s->v0.y, v->y))
243 {
244 if (v->x < s->v0.x)
245 area = 1.0;
246 else
247 area = -1.0;
248 }
249 else
250 area = CROSS(s->v0, s->v1, (*v));
251 }
252 else /* v0 > v1 */
253 {
254 if (FP_EQUAL(s->v1.y, v->y))
255 {
256 if (v->x < s->v1.x)
257 area = 1.0;
258 else
259 area = -1.0;
260 }
261 else if (FP_EQUAL(s->v0.y, v->y))
262 {
263 if (v->x < s->v0.x)
264 area = 1.0;
265 else
266 area = -1.0;
267 }
268 else
269 area = CROSS(s->v1, s->v0, (*v));
270 }
271
272 if (area > 0.0)
273 return TRUE;
274 else
275 return FALSE;
276 }
277
278
279
280 /* Returns true if the corresponding endpoint of the given segment is */
281 /* already inserted into the segment tree. Use the simple test of */
282 /* whether the segment which shares this endpoint is already inserted */
283
284 static int inserted(segnum, whichpt)
285 int segnum;
286 int whichpt;
287 {
288 if (whichpt == FIRSTPT)
289 return seg[seg[segnum].prev].is_inserted;
290 else
291 return seg[seg[segnum].next].is_inserted;
292 }
293
294 /* This is query routine which determines which trapezoid does the
295 * point v lie in. The return value is the trapezoid number.
296 */
297
298 int locate_endpoint(v, vo, r)
299 point_t *v;
300 point_t *vo;
301 int r;
302 {
303 node_t *rptr = &qs[r];
304
305 switch (rptr->nodetype)
306 {
307 case T_SINK:
308 return rptr->trnum;
309
310 case T_Y:
311 if (_greater_than(v, &rptr->yval)) /* above */
312 return locate_endpoint(v, vo, rptr->right);
313 else if (_equal_to(v, &rptr->yval)) /* the point is already */
314 { /* inserted. */
315 if (_greater_than(vo, &rptr->yval)) /* above */
316 return locate_endpoint(v, vo, rptr->right);
317 else
318 return locate_endpoint(v, vo, rptr->left); /* below */
319 }
320 else
321 return locate_endpoint(v, vo, rptr->left); /* below */
322
323 case T_X:
324 if (_equal_to(v, &seg[rptr->segnum].v0) ||
325 _equal_to(v, &seg[rptr->segnum].v1))
326 {
327 if (FP_EQUAL(v->y, vo->y)) /* horizontal segment */
328 {
329 if (vo->x < v->x)
330 return locate_endpoint(v, vo, rptr->left); /* left */
331 else
332 return locate_endpoint(v, vo, rptr->right); /* right */
333 }
334
335 else if (is_left_of(rptr->segnum, vo))
336 return locate_endpoint(v, vo, rptr->left); /* left */
337 else
338 return locate_endpoint(v, vo, rptr->right); /* right */
339 }
340 else if (is_left_of(rptr->segnum, v))
341 return locate_endpoint(v, vo, rptr->left); /* left */
342 else
343 return locate_endpoint(v, vo, rptr->right); /* right */
344
345 default:
346 fprintf(stderr, "Haggu !!!!!\n");
347 break;
348 }
349 }
350
351
352 /* Thread in the segment into the existing trapezoidation. The
353 * limiting trapezoids are given by tfirst and tlast (which are the
354 * trapezoids containing the two endpoints of the segment. Merges all
355 * possible trapezoids which flank this segment and have been recently
356 * divided because of its insertion
357 */
358
359 static int merge_trapezoids(segnum, tfirst, tlast, side)
360 int segnum;
361 int tfirst;
362 int tlast;
363 int side;
364 {
365 int t, tnext, cond;
366 int ptnext;
367
368 /* First merge polys on the LHS */
369 t = tfirst;
370 while ((t > 0) && _greater_than_equal_to(&tr[t].lo, &tr[tlast].lo))
371 {
372 if (side == S_LEFT)
373 cond = ((((tnext = tr[t].d0) > 0) && (tr[tnext].rseg == segnum)) ||
374 (((tnext = tr[t].d1) > 0) && (tr[tnext].rseg == segnum)));
375 else
376 cond = ((((tnext = tr[t].d0) > 0) && (tr[tnext].lseg == segnum)) ||
377 (((tnext = tr[t].d1) > 0) && (tr[tnext].lseg == segnum)));
378
379 if (cond)
380 {
381 if ((tr[t].lseg == tr[tnext].lseg) &&
382 (tr[t].rseg == tr[tnext].rseg)) /* good neighbours */
383 { /* merge them */
384 /* Use the upper node as the new node i.e. t */
385
386 ptnext = qs[tr[tnext].sink].parent;
387
388 if (qs[ptnext].left == tr[tnext].sink)
389 qs[ptnext].left = tr[t].sink;
390 else
391 qs[ptnext].right = tr[t].sink; /* redirect parent */
392
393
394 /* Change the upper neighbours of the lower trapezoids */
395
396 if ((tr[t].d0 = tr[tnext].d0) > 0)
397 if (tr[tr[t].d0].u0 == tnext)
398 tr[tr[t].d0].u0 = t;
399 else if (tr[tr[t].d0].u1 == tnext)
400 tr[tr[t].d0].u1 = t;
401
402 if ((tr[t].d1 = tr[tnext].d1) > 0)
403 if (tr[tr[t].d1].u0 == tnext)
404 tr[tr[t].d1].u0 = t;
405 else if (tr[tr[t].d1].u1 == tnext)
406 tr[tr[t].d1].u1 = t;
407
408 tr[t].lo = tr[tnext].lo;
409 tr[tnext].state = ST_INVALID; /* invalidate the lower */
410 /* trapezium */
411 }
412 else /* not good neighbours */
413 t = tnext;
414 }
415 else /* do not satisfy the outer if */
416 t = tnext;
417
418 } /* end-while */
419
420 return 0;
421 }
422
423
424 /* Add in the new segment into the trapezoidation and update Q and T
425 * structures. First locate the two endpoints of the segment in the
426 * Q-structure. Then start from the topmost trapezoid and go down to
427 * the lower trapezoid dividing all the trapezoids in between .
428 */
429
430 static int add_segment(segnum)
431 int segnum;
432 {
433 segment_t s;
434 segment_t *so = &seg[segnum];
435 int tu, tl, sk, tfirst, tlast, tnext;
436 int tfirstr, tlastr, tfirstl, tlastl;
437 int i1, i2, t, t1, t2, tn;
438 point_t tpt;
439 int tritop = 0, tribot = 0, is_swapped = 0;
440 int tmptriseg;
441
442 s = seg[segnum];
443 if (_greater_than(&s.v1, &s.v0)) /* Get higher vertex in v0 */
444 {
445 int tmp;
446 tpt = s.v0;
447 s.v0 = s.v1;
448 s.v1 = tpt;
449 tmp = s.root0;
450 s.root0 = s.root1;
451 s.root1 = tmp;
452 is_swapped = TRUE;
453 }
454
455 if ((is_swapped) ? !inserted(segnum, LASTPT) :
456 !inserted(segnum, FIRSTPT)) /* insert v0 in the tree */
457 {
458 int tmp_d;
459
460 tu = locate_endpoint(&s.v0, &s.v1, s.root0);
461 tl = newtrap(); /* tl is the new lower trapezoid */
462 tr[tl].state = ST_VALID;
463 tr[tl] = tr[tu];
464 tr[tu].lo.y = tr[tl].hi.y = s.v0.y;
465 tr[tu].lo.x = tr[tl].hi.x = s.v0.x;
466 tr[tu].d0 = tl;
467 tr[tu].d1 = 0;
468 tr[tl].u0 = tu;
469 tr[tl].u1 = 0;
470
471 if (((tmp_d = tr[tl].d0) > 0) && (tr[tmp_d].u0 == tu))
472 tr[tmp_d].u0 = tl;
473 if (((tmp_d = tr[tl].d0) > 0) && (tr[tmp_d].u1 == tu))
474 tr[tmp_d].u1 = tl;
475
476 if (((tmp_d = tr[tl].d1) > 0) && (tr[tmp_d].u0 == tu))
477 tr[tmp_d].u0 = tl;
478 if (((tmp_d = tr[tl].d1) > 0) && (tr[tmp_d].u1 == tu))
479 tr[tmp_d].u1 = tl;
480
481 /* Now update the query structure and obtain the sinks for the */
482 /* two trapezoids */
483
484 i1 = newnode(); /* Upper trapezoid sink */
485 i2 = newnode(); /* Lower trapezoid sink */
486 sk = tr[tu].sink;
487
488 qs[sk].nodetype = T_Y;
489 qs[sk].yval = s.v0;
490 qs[sk].segnum = segnum; /* not really reqd ... maybe later */
491 qs[sk].left = i2;
492 qs[sk].right = i1;
493
494 qs[i1].nodetype = T_SINK;
495 qs[i1].trnum = tu;
496 qs[i1].parent = sk;
497
498 qs[i2].nodetype = T_SINK;
499 qs[i2].trnum = tl;
500 qs[i2].parent = sk;
501
502 tr[tu].sink = i1;
503 tr[tl].sink = i2;
504 tfirst = tl;
505 }
506 else /* v0 already present */
507 { /* Get the topmost intersecting trapezoid */
508 tfirst = locate_endpoint(&s.v0, &s.v1, s.root0);
509 tritop = 1;
510 }
511
512
513 if ((is_swapped) ? !inserted(segnum, FIRSTPT) :
514 !inserted(segnum, LASTPT)) /* insert v1 in the tree */
515 {
516 int tmp_d;
517
518 tu = locate_endpoint(&s.v1, &s.v0, s.root1);
519
520 tl = newtrap(); /* tl is the new lower trapezoid */
521 tr[tl].state = ST_VALID;
522 tr[tl] = tr[tu];
523 tr[tu].lo.y = tr[tl].hi.y = s.v1.y;
524 tr[tu].lo.x = tr[tl].hi.x = s.v1.x;
525 tr[tu].d0 = tl;
526 tr[tu].d1 = 0;
527 tr[tl].u0 = tu;
528 tr[tl].u1 = 0;
529
530 if (((tmp_d = tr[tl].d0) > 0) && (tr[tmp_d].u0 == tu))
531 tr[tmp_d].u0 = tl;
532 if (((tmp_d = tr[tl].d0) > 0) && (tr[tmp_d].u1 == tu))
533 tr[tmp_d].u1 = tl;
534
535 if (((tmp_d = tr[tl].d1) > 0) && (tr[tmp_d].u0 == tu))
536 tr[tmp_d].u0 = tl;
537 if (((tmp_d = tr[tl].d1) > 0) && (tr[tmp_d].u1 == tu))
538 tr[tmp_d].u1 = tl;
539
540 /* Now update the query structure and obtain the sinks for the */
541 /* two trapezoids */
542
543 i1 = newnode(); /* Upper trapezoid sink */
544 i2 = newnode(); /* Lower trapezoid sink */
545 sk = tr[tu].sink;
546
547 qs[sk].nodetype = T_Y;
548 qs[sk].yval = s.v1;
549 qs[sk].segnum = segnum; /* not really reqd ... maybe later */
550 qs[sk].left = i2;
551 qs[sk].right = i1;
552
553 qs[i1].nodetype = T_SINK;
554 qs[i1].trnum = tu;
555 qs[i1].parent = sk;
556
557 qs[i2].nodetype = T_SINK;
558 qs[i2].trnum = tl;
559 qs[i2].parent = sk;
560
561 tr[tu].sink = i1;
562 tr[tl].sink = i2;
563 tlast = tu;
564 }
565 else /* v1 already present */
566 { /* Get the lowermost intersecting trapezoid */
567 tlast = locate_endpoint(&s.v1, &s.v0, s.root1);
568 tribot = 1;
569 }
570
571 /* Thread the segment into the query tree creating a new X-node */
572 /* First, split all the trapezoids which are intersected by s into */
573 /* two */
574
575 t = tfirst; /* topmost trapezoid */
576
577 while ((t > 0) &&
578 _greater_than_equal_to(&tr[t].lo, &tr[tlast].lo))
579 /* traverse from top to bot */
580 {
581 int t_sav, tn_sav;
582 sk = tr[t].sink;
583 i1 = newnode(); /* left trapezoid sink */
584 i2 = newnode(); /* right trapezoid sink */
585
586 qs[sk].nodetype = T_X;
587 qs[sk].segnum = segnum;
588 qs[sk].left = i1;
589 qs[sk].right = i2;
590
591 qs[i1].nodetype = T_SINK; /* left trapezoid (use existing one) */
592 qs[i1].trnum = t;
593 qs[i1].parent = sk;
594
595 qs[i2].nodetype = T_SINK; /* right trapezoid (allocate new) */
596 qs[i2].trnum = tn = newtrap();
597 tr[tn].state = ST_VALID;
598 qs[i2].parent = sk;
599
600 if (t == tfirst)
601 tfirstr = tn;
602 if (_equal_to(&tr[t].lo, &tr[tlast].lo))
603 tlastr = tn;
604
605 tr[tn] = tr[t];
606 tr[t].sink = i1;
607 tr[tn].sink = i2;
608 t_sav = t;
609 tn_sav = tn;
610
611 /* error */
612
613 if ((tr[t].d0 <= 0) && (tr[t].d1 <= 0)) /* case cannot arise */
614 {
615 fprintf(stderr, "add_segment: error\n");
616 break;
617 }
618
619 /* only one trapezoid below. partition t into two and make the */
620 /* two resulting trapezoids t and tn as the upper neighbours of */
621 /* the sole lower trapezoid */
622
623 else if ((tr[t].d0 > 0) && (tr[t].d1 <= 0))
624 { /* Only one trapezoid below */
625 if ((tr[t].u0 > 0) && (tr[t].u1 > 0))
626 { /* continuation of a chain from abv. */
627 if (tr[t].usave > 0) /* three upper neighbours */
628 {
629 if (tr[t].uside == S_LEFT)
630 {
631 tr[tn].u0 = tr[t].u1;
632 tr[t].u1 = -1;
633 tr[tn].u1 = tr[t].usave;
634
635 tr[tr[t].u0].d0 = t;
636 tr[tr[tn].u0].d0 = tn;
637 tr[tr[tn].u1].d0 = tn;
638 }
639 else /* intersects in the right */
640 {
641 tr[tn].u1 = -1;
642 tr[tn].u0 = tr[t].u1;
643 tr[t].u1 = tr[t].u0;
644 tr[t].u0 = tr[t].usave;
645
646 tr[tr[t].u0].d0 = t;
647 tr[tr[t].u1].d0 = t;
648 tr[tr[tn].u0].d0 = tn;
649 }
650
651 tr[t].usave = tr[tn].usave = 0;
652 }
653 else /* No usave.... simple case */
654 {
655 tr[tn].u0 = tr[t].u1;
656 tr[t].u1 = tr[tn].u1 = -1;
657 tr[tr[tn].u0].d0 = tn;
658 }
659 }
660 else
661 { /* fresh seg. or upward cusp */
662 int tmp_u = tr[t].u0;
663 int td0, td1;
664 if (((td0 = tr[tmp_u].d0) > 0) &&
665 ((td1 = tr[tmp_u].d1) > 0))
666 { /* upward cusp */
667 if ((tr[td0].rseg > 0) &&
668 !is_left_of(tr[td0].rseg, &s.v1))
669 {
670 tr[t].u0 = tr[t].u1 = tr[tn].u1 = -1;
671 tr[tr[tn].u0].d1 = tn;
672 }
673 else /* cusp going leftwards */
674 {
675 tr[tn].u0 = tr[tn].u1 = tr[t].u1 = -1;
676 tr[tr[t].u0].d0 = t;
677 }
678 }
679 else /* fresh segment */
680 {
681 tr[tr[t].u0].d0 = t;
682 tr[tr[t].u0].d1 = tn;
683 }
684 }
685
686 if (FP_EQUAL(tr[t].lo.y, tr[tlast].lo.y) &&
687 FP_EQUAL(tr[t].lo.x, tr[tlast].lo.x) && tribot)
688 { /* bottom forms a triangle */
689
690 if (is_swapped)
691 tmptriseg = seg[segnum].prev;
692 else
693 tmptriseg = seg[segnum].next;
694
695 if ((tmptriseg > 0) && is_left_of(tmptriseg, &s.v0))
696 {
697 /* L-R downward cusp */
698 tr[tr[t].d0].u0 = t;
699 tr[tn].d0 = tr[tn].d1 = -1;
700 }
701 else
702 {
703 /* R-L downward cusp */
704 tr[tr[tn].d0].u1 = tn;
705 tr[t].d0 = tr[t].d1 = -1;
706 }
707 }
708 else
709 {
710 if ((tr[tr[t].d0].u0 > 0) && (tr[tr[t].d0].u1 > 0))
711 {
712 if (tr[tr[t].d0].u0 == t) /* passes thru LHS */
713 {
714 tr[tr[t].d0].usave = tr[tr[t].d0].u1;
715 tr[tr[t].d0].uside = S_LEFT;
716 }
717 else
718 {
719 tr[tr[t].d0].usave = tr[tr[t].d0].u0;
720 tr[tr[t].d0].uside = S_RIGHT;
721 }
722 }
723 tr[tr[t].d0].u0 = t;
724 tr[tr[t].d0].u1 = tn;
725 }
726
727 t = tr[t].d0;
728 }
729
730
731 else if ((tr[t].d0 <= 0) && (tr[t].d1 > 0))
732 { /* Only one trapezoid below */
733 if ((tr[t].u0 > 0) && (tr[t].u1 > 0))
734 { /* continuation of a chain from abv. */
735 if (tr[t].usave > 0) /* three upper neighbours */
736 {
737 if (tr[t].uside == S_LEFT)
738 {
739 tr[tn].u0 = tr[t].u1;
740 tr[t].u1 = -1;
741 tr[tn].u1 = tr[t].usave;
742
743 tr[tr[t].u0].d0 = t;
744 tr[tr[tn].u0].d0 = tn;
745 tr[tr[tn].u1].d0 = tn;
746 }
747 else /* intersects in the right */
748 {
749 tr[tn].u1 = -1;
750 tr[tn].u0 = tr[t].u1;
751 tr[t].u1 = tr[t].u0;
752 tr[t].u0 = tr[t].usave;
753
754 tr[tr[t].u0].d0 = t;
755 tr[tr[t].u1].d0 = t;
756 tr[tr[tn].u0].d0 = tn;
757 }
758
759 tr[t].usave = tr[tn].usave = 0;
760 }
761 else /* No usave.... simple case */
762 {
763 tr[tn].u0 = tr[t].u1;
764 tr[t].u1 = tr[tn].u1 = -1;
765 tr[tr[tn].u0].d0 = tn;
766 }
767 }
768 else
769 { /* fresh seg. or upward cusp */
770 int tmp_u = tr[t].u0;
771 int td0, td1;
772 if (((td0 = tr[tmp_u].d0) > 0) &&
773 ((td1 = tr[tmp_u].d1) > 0))
774 { /* upward cusp */
775 if ((tr[td0].rseg > 0) &&
776 !is_left_of(tr[td0].rseg, &s.v1))
777 {
778 tr[t].u0 = tr[t].u1 = tr[tn].u1 = -1;
779 tr[tr[tn].u0].d1 = tn;
780 }
781 else
782 {
783 tr[tn].u0 = tr[tn].u1 = tr[t].u1 = -1;
784 tr[tr[t].u0].d0 = t;
785 }
786 }
787 else /* fresh segment */
788 {
789 tr[tr[t].u0].d0 = t;
790 tr[tr[t].u0].d1 = tn;
791 }
792 }
793
794 if (FP_EQUAL(tr[t].lo.y, tr[tlast].lo.y) &&
795 FP_EQUAL(tr[t].lo.x, tr[tlast].lo.x) && tribot)
796 { /* bottom forms a triangle */
797 int tmpseg;
798
799 if (is_swapped)
800 tmptriseg = seg[segnum].prev;
801 else
802 tmptriseg = seg[segnum].next;
803
804 if ((tmpseg > 0) && is_left_of(tmpseg, &s.v0))
805 {
806 /* L-R downward cusp */
807 tr[tr[t].d1].u0 = t;
808 tr[tn].d0 = tr[tn].d1 = -1;
809 }
810 else
811 {
812 /* R-L downward cusp */
813 tr[tr[tn].d1].u1 = tn;
814 tr[t].d0 = tr[t].d1 = -1;
815 }
816 }
817 else
818 {
819 if ((tr[tr[t].d1].u0 > 0) && (tr[tr[t].d1].u1 > 0))
820 {
821 if (tr[tr[t].d1].u0 == t) /* passes thru LHS */
822 {
823 tr[tr[t].d1].usave = tr[tr[t].d1].u1;
824 tr[tr[t].d1].uside = S_LEFT;
825 }
826 else
827 {
828 tr[tr[t].d1].usave = tr[tr[t].d1].u0;
829 tr[tr[t].d1].uside = S_RIGHT;
830 }
831 }
832 tr[tr[t].d1].u0 = t;
833 tr[tr[t].d1].u1 = tn;
834 }
835
836 t = tr[t].d1;
837 }
838
839 /* two trapezoids below. Find out which one is intersected by */
840 /* this segment and proceed down that one */
841
842 else
843 {
844 int tmpseg = tr[tr[t].d0].rseg;
845 double y0, yt;
846 point_t tmppt;
847 int tnext, i_d0, i_d1;
848
849 i_d0 = i_d1 = FALSE;
850 if (FP_EQUAL(tr[t].lo.y, s.v0.y))
851 {
852 if (tr[t].lo.x > s.v0.x)
853 i_d0 = TRUE;
854 else
855 i_d1 = TRUE;
856 }
857 else
858 {
859 tmppt.y = y0 = tr[t].lo.y;
860 yt = (y0 - s.v0.y)/(s.v1.y - s.v0.y);
861 tmppt.x = s.v0.x + yt * (s.v1.x - s.v0.x);
862
863 if (_less_than(&tmppt, &tr[t].lo))
864 i_d0 = TRUE;
865 else
866 i_d1 = TRUE;
867 }
868
869 /* check continuity from the top so that the lower-neighbour */
870 /* values are properly filled for the upper trapezoid */
871
872 if ((tr[t].u0 > 0) && (tr[t].u1 > 0))
873 { /* continuation of a chain from abv. */
874 if (tr[t].usave > 0) /* three upper neighbours */
875 {
876 if (tr[t].uside == S_LEFT)
877 {
878 tr[tn].u0 = tr[t].u1;
879 tr[t].u1 = -1;
880 tr[tn].u1 = tr[t].usave;
881
882 tr[tr[t].u0].d0 = t;
883 tr[tr[tn].u0].d0 = tn;
884 tr[tr[tn].u1].d0 = tn;
885 }
886 else /* intersects in the right */
887 {
888 tr[tn].u1 = -1;
889 tr[tn].u0 = tr[t].u1;
890 tr[t].u1 = tr[t].u0;
891 tr[t].u0 = tr[t].usave;
892
893 tr[tr[t].u0].d0 = t;
894 tr[tr[t].u1].d0 = t;
895 tr[tr[tn].u0].d0 = tn;
896 }
897
898 tr[t].usave = tr[tn].usave = 0;
899 }
900 else /* No usave.... simple case */
901 {
902 tr[tn].u0 = tr[t].u1;
903 tr[tn].u1 = -1;
904 tr[t].u1 = -1;
905 tr[tr[tn].u0].d0 = tn;
906 }
907 }
908 else
909 { /* fresh seg. or upward cusp */
910 int tmp_u = tr[t].u0;
911 int td0, td1;
912 if (((td0 = tr[tmp_u].d0) > 0) &&
913 ((td1 = tr[tmp_u].d1) > 0))
914 { /* upward cusp */
915 if ((tr[td0].rseg > 0) &&
916 !is_left_of(tr[td0].rseg, &s.v1))
917 {
918 tr[t].u0 = tr[t].u1 = tr[tn].u1 = -1;
919 tr[tr[tn].u0].d1 = tn;
920 }
921 else
922 {
923 tr[tn].u0 = tr[tn].u1 = tr[t].u1 = -1;
924 tr[tr[t].u0].d0 = t;
925 }
926 }
927 else /* fresh segment */
928 {
929 tr[tr[t].u0].d0 = t;
930 tr[tr[t].u0].d1 = tn;
931 }
932 }
933
934 if (FP_EQUAL(tr[t].lo.y, tr[tlast].lo.y) &&
935 FP_EQUAL(tr[t].lo.x, tr[tlast].lo.x) && tribot)
936 {
937 /* this case arises only at the lowest trapezoid.. i.e.
938 tlast, if the lower endpoint of the segment is
939 already inserted in the structure */
940
941 tr[tr[t].d0].u0 = t;
942 tr[tr[t].d0].u1 = -1;
943 tr[tr[t].d1].u0 = tn;
944 tr[tr[t].d1].u1 = -1;
945
946 tr[tn].d0 = tr[t].d1;
947 tr[t].d1 = tr[tn].d1 = -1;
948
949 tnext = tr[t].d1;
950 }
951 else if (i_d0)
952 /* intersecting d0 */
953 {
954 tr[tr[t].d0].u0 = t;
955 tr[tr[t].d0].u1 = tn;
956 tr[tr[t].d1].u0 = tn;
957 tr[tr[t].d1].u1 = -1;
958
959 /* new code to determine the bottom neighbours of the */
960 /* newly partitioned trapezoid */
961
962 tr[t].d1 = -1;
963
964 tnext = tr[t].d0;
965 }
966 else /* intersecting d1 */
967 {
968 tr[tr[t].d0].u0 = t;
969 tr[tr[t].d0].u1 = -1;
970 tr[tr[t].d1].u0 = t;
971 tr[tr[t].d1].u1 = tn;
972
973 /* new code to determine the bottom neighbours of the */
974 /* newly partitioned trapezoid */
975
976 tr[tn].d0 = tr[t].d1;
977 tr[tn].d1 = -1;
978
979 tnext = tr[t].d1;
980 }
981
982 t = tnext;
983 }
984
985 tr[t_sav].rseg = tr[tn_sav].lseg = segnum;
986 } /* end-while */
987
988 /* Now combine those trapezoids which share common segments. We can */
989 /* use the pointers to the parent to connect these together. This */
990 /* works only because all these new trapezoids have been formed */
991 /* due to splitting by the segment, and hence have only one parent */
992
993 tfirstl = tfirst;
994 tlastl = tlast;
995 merge_trapezoids(segnum, tfirstl, tlastl, S_LEFT);
996 merge_trapezoids(segnum, tfirstr, tlastr, S_RIGHT);
997
998 seg[segnum].is_inserted = TRUE;
999 return 0;
1000 }
1001
1002
1003 /* Update the roots stored for each of the endpoints of the segment.
1004 * This is done to speed up the location-query for the endpoint when
1005 * the segment is inserted into the trapezoidation subsequently
1006 */
1007 static int find_new_roots(segnum)
1008 int segnum;
1009 {
1010 segment_t *s = &seg[segnum];
1011
1012 if (s->is_inserted)
1013 return 0;
1014
1015 s->root0 = locate_endpoint(&s->v0, &s->v1, s->root0);
1016 s->root0 = tr[s->root0].sink;
1017
1018 s->root1 = locate_endpoint(&s->v1, &s->v0, s->root1);
1019 s->root1 = tr[s->root1].sink;
1020 return 0;
1021 }
1022
1023
1024 /* Main routine to perform trapezoidation */
1025 int construct_trapezoids(nseg)
1026 int nseg;
1027 {
1028 register int i;
1029 int root, h;
1030
1031 /* Add the first segment and get the query structure and trapezoid */
1032 /* list initialised */
1033
1034 root = init_query_structure(choose_segment());
1035
1036 for (i = 1; i <= nseg; i++)
1037 seg[i].root0 = seg[i].root1 = root;
1038
1039 for (h = 1; h <= math_logstar_n(nseg); h++)
1040 {
1041 for (i = math_N(nseg, h -1) + 1; i <= math_N(nseg, h); i++)
1042 add_segment(choose_segment());
1043
1044 /* Find a new root for each of the segment endpoints */
1045 for (i = 1; i <= nseg; i++)
1046 find_new_roots(i);
1047 }
1048
1049 for (i = math_N(nseg, math_logstar_n(nseg)) + 1; i <= nseg; i++)
1050 add_segment(choose_segment());
1051
1052 return 0;
1053 }
1054
1055

   
Visit the ZANavi Wiki