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 |
|