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