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