src/polygon1.c

   1 /*
   2        polygon clipping functions. harry eaton implemented the algorithm
   3        described in "A Closed Set of Algorithms for Performing Set
   4        Operations on Polygonal Regions in the Plane" which the original
   5        code did not do. I also modified it for integer coordinates
   6        and faster computation. The license for this modified copy was
   7        switched to the GPL per term (3) of the original LGPL license.
   8        Copyright (C) 2006 harry eaton
   9
  10    based on:
  11        poly_Boolean: a polygon clip library
  12        Copyright (C) 1997  Alexey Nikitin, Michael Leonov
  13        (also the authors of the paper describing the actual algorithm)
  14        leonov@propro.iis.nsk.su
  15
  16    in turn based on:
  17        nclip: a polygon clip library
  18        Copyright (C) 1993  Klamer Schutte
  19
  20        This program is free software; you can redistribute it and/or
  21        modify it under the terms of the GNU General Public
  22        License as published by the Free Software Foundation; either
  23        version 2 of the License, or (at your option) any later version.
  24
  25        This program is distributed in the hope that it will be useful,
  26        but WITHOUT ANY WARRANTY; without even the implied warranty of
  27        MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  28        General Public License for more details.
  29
  30        You should have received a copy of the GNU General Public
  31        License along with this program; if not, write to the Free
  32        Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  33
  34       polygon1.c
  35       (C) 1997 Alexey Nikitin, Michael Leonov
  36       (C) 1993 Klamer Schutte
  37
  38       all cases where original (Klamer Schutte) code is present
  39       are marked
  40 */
  41
  42 #include        <assert.h>
  43 #include        <stdlib.h>
  44 #include        <stdio.h>
  45 #include        <setjmp.h>
  46 #include        <math.h>
  47 #include        <string.h>
  48
  49 #include "global.h"
  50 #include "rtree.h"
  51 #include "heap.h"
  52
  53 #define ROUND(a) (long)((a) > 0 ? ((a) + 0.5) : ((a) - 0.5))
  54
  55 #define EPSILON (1E-8)
  56 #define IsZero(a, b) (fabs((a) - (b)) < EPSILON)
  57 #ifndef ABS
  58 #define ABS(x) ((x) < 0 ? -(x) : (x))
  59 #endif
  60
  61 /*********************************************************************/
  62 /*              L o n g   V e c t o r   S t u f f                    */
  63 /*********************************************************************/
  64
  65 #define Vcopy(a,b) {(a)[0]=(b)[0];(a)[1]=(b)[1];}
  66 int vect_equal (Vector v1, Vector v2);
  67 void vect_init (Vector v, double x, double y);
  68 void vect_sub (Vector res, Vector v2, Vector v3);
  69
  70 void vect_min (Vector res, Vector v2, Vector v3);
  71 void vect_max (Vector res, Vector v2, Vector v3);
  72
  73 double vect_dist2 (Vector v1, Vector v2);
  74 double vect_det2 (Vector v1, Vector v2);
  75 double vect_len2 (Vector v1);
  76
  77 int vect_inters2 (Vector A, Vector B, Vector C, Vector D, Vector S1,
  78                   Vector S2);
  79
  80 /* note that a vertex v's Flags.status represents the edge defined by
  81  * v to v->next (i.e. the edge is forward of v)
  82  */
  83 #define ISECTED 3
  84 #define UNKNWN  0
  85 #define INSIDE  1
  86 #define OUTSIDE 2
  87 #define SHARED 3
  88 #define SHARED2 4
  89
  90 #define TOUCHES 99
  91
  92 #define NODE_LABEL(n)  ((n)->Flags.status)
  93 #define LABEL_NODE(n,l) ((n)->Flags.status = (l))
  94
  95 #define error(code)  longjmp(*(e), code)
  96
  97 #define MemGet(ptr, type) \
  98 if (UNLIKELY (((ptr) = malloc(sizeof(type))) == NULL)) \
  99     error(err_no_memory);
 100
 101 #undef DEBUG_LABEL
 102 #undef DEBUG_ALL_LABELS
 103 #undef DEBUG_JUMP
 104 #undef DEBUG_GATHER
 105 #undef DEBUG_ANGLE
 106 #undef DEBUG
 107 #ifdef DEBUG
 108 #define DEBUGP(...) fprintf(stderr, ## __VA_ARGS__)
 109 #else
 110 #define DEBUGP(...)
 111 #endif
 112
 113 /* ///////////////////////////////////////////////////////////////////////////// * /
 114 / *  2-Dimentional stuff
 115 / * ///////////////////////////////////////////////////////////////////////////// */
 116
 117 #define Vsub2(r,a,b)    {(r)[0] = (a)[0] - (b)[0]; (r)[1] = (a)[1] - (b)[1];}
 118 #define Vadd2(r,a,b)    {(r)[0] = (a)[0] + (b)[0]; (r)[1] = (a)[1] + (b)[1];}
 119 #define Vsca2(r,a,s)    {(r)[0] = (a)[0] * (s); (r)[1] = (a)[1] * (s);}
 120 #define Vcpy2(r,a)              {(r)[0] = (a)[0]; (r)[1] = (a)[1];}
 121 #define Vequ2(a,b)              ((a)[0] == (b)[0] && (a)[1] == (b)[1])
 122 #define Vadds(r,a,b,s)  {(r)[0] = ((a)[0] + (b)[0]) * (s); (r)[1] = ((a)[1] + (b)[1]) * (s);}
 123 #define Vswp2(a,b) { long t; \
 124         t = (a)[0], (a)[0] = (b)[0], (b)[0] = t; \
 125         t = (a)[1], (a)[1] = (b)[1], (b)[1] = t; \
 126 }
 127
 128 #ifdef DEBUG
 129 static char *theState (VNODE * v);
 130
 131 static void
 132 pline_dump (VNODE * v)
 133 {
 134   VNODE *s;
 135
 136   s = v;
 137   do
 138     {
 139       if (v != s)
 140         fprintf (stderr, "%d %d 10 10 \"%s\"]\n", v->point[0], v->point[1],
 141                  theState (v));
 142       fprintf (stderr, "Line [%d %d ", v->point[0], v->point[1]);
 143     }
 144   while ((v = v->next) != s);
 145   fprintf (stderr, "%d %d 10 10 \"%s\"]\n", v->point[0], v->point[1],
 146            theState (v));
 147   fprintf (stderr, "NEXT PLINE\n");
 148 }
 149
 150 static void
 151 poly_dump (POLYAREA * p)
 152 {
 153   POLYAREA *f = p;
 154
 155   do
 156     {
 157       pline_dump (&p->contours->head);
 158     }
 159   while ((p = p->f) != f);
 160   fprintf (stderr, "NEXT_POLY\n");
 161 }
 162 #endif
 163
 164 /***************************************************************/
 165 /* routines for processing intersections */
 166
 167 /*
 168 node_add
 169  (C) 1993 Klamer Schutte
 170  (C) 1997 Alexey Nikitin, Michael Leonov
 171  (C) 2006 harry eaton
 172
 173  returns a bit field in new_point that indicates where the
 174  point was.
 175  1 means a new node was created and inserted
 176  4 means the intersection was not on the dest point
 177 */
 178 static VNODE *
 179 node_add (VNODE * dest, Vector po, int *new_point)
 180 {
 181   VNODE *p;
 182
 183   if (vect_equal (po, dest->point))
 184     return dest;
 185   if (vect_equal (po, dest->next->point))
 186     {
 187       (*new_point) += 4;
 188       return dest->next;
 189     }
 190   p = poly_CreateNode (po);
 191   if (p == NULL)
 192     return NULL;
 193   (*new_point) += 5;
 194   p->prev = dest;
 195   p->next = dest->next;
 196   p->cvc_prev = p->cvc_next = NULL;
 197   p->Flags.status = UNKNWN;
 198   return (dest->next = dest->next->prev = p);
 199 }                               /* node_add */
 200
 201 #define ISECT_BAD_PARAM (-1)
 202 #define ISECT_NO_MEMORY (-2)
 203
 204
 205 /*
 206 new_descriptor
 207   (C) 2006 harry eaton
 208 */
 209 static CVCList *
 210 new_descriptor (VNODE * a, char poly, char side)
 211 {
 212   CVCList *l = (CVCList *) malloc (sizeof (CVCList));
 213   Vector v;
 214   register double ang, dx, dy;
 215
 216   if (!l)
 217     return NULL;
 218   l->head = NULL;
 219   l->parent = a;
 220   l->poly = poly;
 221   l->side = side;
 222   l->next = l->prev = l;
 223   if (side == 'P')              /* previous */
 224     vect_sub (v, a->prev->point, a->point);
 225   else                          /* next */
 226     vect_sub (v, a->next->point, a->point);
 227   /* Uses slope/(slope+1) in quadrant 1 as a proxy for the angle.
 228    * It still has the same monotonic sort result
 229    * and is far less expensive to compute than the real angle.
 230    */
 231   if (vect_equal (v, vect_zero))
 232     {
 233       if (side == 'P')
 234         {
 235           if (a->prev->cvc_prev == (CVCList *) - 1)
 236             a->prev->cvc_prev = a->prev->cvc_next = NULL;
 237           poly_ExclVertex (a->prev);
 238           vect_sub (v, a->prev->point, a->point);
 239         }
 240       else
 241         {
 242           if (a->next->cvc_prev == (CVCList *) - 1)
 243             a->next->cvc_prev = a->next->cvc_next = NULL;
 244           poly_ExclVertex (a->next);
 245           vect_sub (v, a->next->point, a->point);
 246         }
 247     }
 248   assert (!vect_equal (v, vect_zero));
 249   dx = fabs ((double) v[0]);
 250   dy = fabs ((double) v[1]);
 251   ang = dy / (dy + dx);
 252   /* now move to the actual quadrant */
 253   if (v[0] < 0 && v[1] >= 0)
 254     ang = 2.0 - ang;            /* 2nd quadrant */
 255   else if (v[0] < 0 && v[1] < 0)
 256     ang += 2.0;                 /* 3rd quadrant */
 257   else if (v[0] >= 0 && v[1] < 0)
 258     ang = 4.0 - ang;            /* 4th quadrant */
 259   l->angle = ang;
 260   assert (ang >= 0.0 && ang <= 4.0);
 261 #ifdef DEBUG_ANGLE
 262   DEBUGP ("node on %c at (%ld,%ld) assigned angle %g on side %c\n", poly,
 263           a->point[0], a->point[1], ang, side);
 264 #endif
 265   return l;
 266 }
 267
 268 /*
 269 insert_descriptor
 270   (C) 2006 harry eaton
 271
 272    argument a is a cross-vertex node.
 273    argument poly is the polygon it comes from ('A' or 'B')
 274    argument side is the side this descriptor goes on ('P' for previous
 275    'N' for next.
 276    argument start is the head of the list of cvclists
 277 */
 278 static CVCList *
 279 insert_descriptor (VNODE * a, char poly, char side, CVCList * start)
 280 {
 281   CVCList *l, *new, *big, *small;
 282
 283   if (!(new = new_descriptor (a, poly, side)))
 284     return NULL;
 285   /* search for the CVCList for this point */
 286   if (!start)
 287     {
 288       start = new;              /* return is also new, so we know where start is */
 289       start->head = new;        /* circular list */
 290       return new;
 291     }
 292   else
 293     {
 294       l = start;
 295       do
 296         {
 297           assert (l->head);
 298           if (l->parent->point[0] == a->point[0]
 299               && l->parent->point[1] == a->point[1])
 300             {                   /* this CVCList is at our point */
 301               start = l;
 302               new->head = l->head;
 303               break;
 304             }
 305           if (l->head->parent->point[0] == start->parent->point[0]
 306               && l->head->parent->point[1] == start->parent->point[1])
 307             {
 308               /* this seems to be a new point */
 309               /* link this cvclist to the list of all cvclists */
 310               for (; l->head != new; l = l->next)
 311                 l->head = new;
 312               new->head = start;
 313               return new;
 314             }
 315           l = l->head;
 316         }
 317       while (1);
 318     }
 319   assert (start);
 320   l = big = small = start;
 321   do
 322     {
 323       if (l->next->angle < l->angle)    /* find start/end of list */
 324         {
 325           small = l->next;
 326           big = l;
 327         }
 328       else if (new->angle >= l->angle && new->angle <= l->next->angle)
 329         {
 330           /* insert new cvc if it lies between existing points */
 331           new->prev = l;
 332           new->next = l->next;
 333           l->next = l->next->prev = new;
 334           return new;
 335         }
 336     }
 337   while ((l = l->next) != start);
 338   /* didn't find it between points, it must go on an end */
 339   if (big->angle <= new->angle)
 340     {
 341       new->prev = big;
 342       new->next = big->next;
 343       big->next = big->next->prev = new;
 344       return new;
 345     }
 346   assert (small->angle >= new->angle);
 347   new->next = small;
 348   new->prev = small->prev;
 349   small->prev = small->prev->next = new;
 350   return new;
 351 }
 352
 353 /*
 354 node_add_point
 355  (C) 1993 Klamer Schutte
 356  (C) 1997 Alexey Nikitin, Michael Leonov
 357
 358  return 1 if new node in b, 2 if new node in a and 3 if new node in both
 359 */
 360
 361 static int
 362 node_add_point (VNODE * a, VNODE * b, Vector p)
 363 {
 364   int res = 0;
 365
 366   VNODE *node_a, *node_b;
 367
 368   node_a = node_add (a, p, &res);
 369   res += res;
 370   node_b = node_add (b, p, &res);
 371
 372   if (node_a == NULL || node_b == NULL)
 373     return ISECT_NO_MEMORY;
 374   node_b->cvc_prev = node_b->cvc_next = (CVCList *) - 1;
 375   node_a->cvc_prev = node_a->cvc_next = (CVCList *) - 1;
 376   return res;
 377 }                               /* node_add_point */
 378
 379 /*
 380 node_label
 381  (C) 2006 harry eaton
 382 */
 383 static unsigned int
 384 node_label (VNODE * pn)
 385 {
 386   CVCList *l;
 387   char this_poly;
 388   int region = UNKNWN;
 389
 390   assert (pn);
 391   assert (pn->cvc_prev);
 392   this_poly = pn->cvc_prev->poly;
 393   /* search clockwise in the cross vertex connectivity (CVC) list
 394    *
 395    * check for shared edges, and check if our edges
 396    * are ever found between the other poly's entry and exit
 397    */
 398 #ifdef DEBUG_LABEL
 399   DEBUGP ("CVCLIST for point (%ld,%ld)\n", pn->point[0], pn->point[1]);
 400 #endif
 401   /* first find whether we're starting inside or outside */
 402   for (l = pn->cvc_prev->prev; l != pn->cvc_prev; l = l->prev)
 403     {
 404       if (l->poly != this_poly)
 405         {
 406           if (l->side == 'P')
 407             region = INSIDE;
 408           else
 409             region = OUTSIDE;
 410         }
 411     }
 412   l = pn->cvc_prev;
 413   do
 414     {
 415       assert (l->angle >= 0 && l->angle <= 4.0);
 416 #ifdef DEBUG_LABEL
 417       DEBUGP ("  poly %c side %c angle = %g\n", l->poly, l->side, l->angle);
 418 #endif
 419       if (l->poly != this_poly)
 420         {
 421           if (l->side == 'P')
 422             {
 423               region = INSIDE;
 424               if (l->parent->prev->point[0] == pn->prev->point[0] &&
 425                   l->parent->prev->point[1] == pn->prev->point[1])
 426                 {
 427                   LABEL_NODE (pn->prev, SHARED);        /* incoming is shared */
 428                   pn->prev->shared = l->parent->prev;
 429                 }
 430               else if (l->parent->prev->point[0] == pn->next->point[0] &&
 431                        l->parent->prev->point[1] == pn->next->point[1])
 432                 {
 433                   LABEL_NODE (pn, SHARED2);     /* outgoing is shared2 */
 434                   pn->shared = l->parent->prev;
 435                 }
 436             }
 437           else
 438             {
 439               region = OUTSIDE;
 440               if (l->parent->next->point[0] == pn->next->point[0] &&
 441                   l->parent->next->point[1] == pn->next->point[1])
 442                 {
 443                   LABEL_NODE (pn, SHARED);
 444                   pn->shared = l->parent;
 445                 }
 446               else if (l->parent->next->point[0] == pn->prev->point[0] &&
 447                        l->parent->next->point[1] == pn->prev->point[1])
 448                 {
 449                   LABEL_NODE (pn->prev, SHARED2);       /* outgoing is shared2 */
 450                   pn->prev->shared = l->parent;
 451                 }
 452             }
 453         }
 454       else
 455         {
 456           VNODE *v;
 457           if (l->side == 'P')
 458             v = l->parent->prev;
 459           else
 460             v = l->parent;
 461           if (NODE_LABEL (v) != SHARED && NODE_LABEL (v) != SHARED2)
 462             {
 463 #ifdef DEBUG_LABEL
 464               /* debugging */
 465               if (NODE_LABEL (v) != UNKNWN && NODE_LABEL (v) != region)
 466                 {
 467                   CVCList *x = l;
 468                   LABEL_NODE (v, region);
 469                   pline_dump (v);
 470                   do
 471                     {
 472                       fprintf (stderr, "poly %c\n", x->poly);
 473                       pline_dump (x->parent);
 474                     }
 475                   while ((x = x->next) != l);
 476                 }
 477 #endif
 478               assert (NODE_LABEL (v) == UNKNWN || NODE_LABEL (v) == region);
 479               LABEL_NODE (v, region);
 480             }
 481         }
 482     }
 483   while ((l = l->prev) != pn->cvc_prev);
 484 #ifdef DEBUG_LABEL
 485   DEBUGP ("\n");
 486 #endif
 487   assert (NODE_LABEL (pn) != UNKNWN && NODE_LABEL (pn->prev) != UNKNWN);
 488   if (NODE_LABEL (pn) == UNKNWN || NODE_LABEL (pn->prev) == UNKNWN)
 489     return UNKNWN;
 490   if (NODE_LABEL (pn) == INSIDE || NODE_LABEL (pn) == OUTSIDE)
 491     return NODE_LABEL (pn);
 492   return UNKNWN;
 493 }                               /* node_label */
 494
 495 /*
 496  add_descriptors
 497  (C) 2006 harry eaton
 498 */
 499 static CVCList *
 500 add_descriptors (PLINE * pl, char poly, CVCList * list)
 501 {
 502   VNODE *node = &pl->head;
 503
 504   do
 505     {
 506       if (node->cvc_prev)
 507         {
 508           assert (node->cvc_prev == (CVCList *) - 1
 509                   && node->cvc_next == (CVCList *) - 1);
 510           list = node->cvc_prev = insert_descriptor (node, poly, 'P', list);
 511           if (!node->cvc_prev)
 512             return NULL;
 513           list = node->cvc_next = insert_descriptor (node, poly, 'N', list);
 514           if (!node->cvc_next)
 515             return NULL;
 516         }
 517     }
 518   while ((node = node->next) != &pl->head);
 519   return list;
 520 }
 521
 522 static inline void
 523 cntrbox_adjust (PLINE * c, Vector p)
 524 {
 525   c->xmin = min (c->xmin, p[0]);
 526   c->xmax = max (c->xmax, p[0] + 1);
 527   c->ymin = min (c->ymin, p[1]);
 528   c->ymax = max (c->ymax, p[1] + 1);
 529 }
 530
 531 /* some structures for handling segment intersections using the rtrees */
 532
 533 typedef struct seg
 534 {
 535   BoxType box;
 536   VNODE *v;
 537   PLINE *p;
 538 } seg;
 539
 540 typedef struct info
 541 {
 542   double m, b;
 543   rtree_t *tree;
 544   VNODE *v;
 545   struct seg *s;
 546   jmp_buf env, sego, *touch;
 547 } info;
 548
 549 /*
 550  * adjust_tree()
 551  * (C) 2006 harry eaton
 552  * This replaces the segment in the tree with the two new segments after
 553  * a vertex has been added
 554  */
 555 static int
 556 adjust_tree (rtree_t * tree, struct seg *s)
 557 {
 558   struct seg *q;
 559
 560   q = malloc (sizeof (struct seg));
 561   if (!q)
 562     return 1;
 563   q->v = s->v;
 564   q->p = s->p;
 565   q->box.X1 = min (q->v->point[0], q->v->next->point[0]);
 566   q->box.X2 = max (q->v->point[0], q->v->next->point[0]) + 1;
 567   q->box.Y1 = min (q->v->point[1], q->v->next->point[1]);
 568   q->box.Y2 = max (q->v->point[1], q->v->next->point[1]) + 1;
 569   r_insert_entry (tree, (const BoxType *) q, 1);
 570   q = malloc (sizeof (struct seg));
 571   if (!q)
 572     return 1;
 573   q->v = s->v->next;
 574   q->p = s->p;
 575   q->box.X1 = min (q->v->point[0], q->v->next->point[0]);
 576   q->box.X2 = max (q->v->point[0], q->v->next->point[0]) + 1;
 577   q->box.Y1 = min (q->v->point[1], q->v->next->point[1]);
 578   q->box.Y2 = max (q->v->point[1], q->v->next->point[1]) + 1;
 579   r_insert_entry (tree, (const BoxType *) q, 1);
 580   r_delete_entry (tree, (const BoxType *) s);
 581   return 0;
 582 }
 583
 584 /*
 585  * seg_in_region()
 586  * (C) 2006, harry eaton
 587  * This prunes the search for boxes that don't intersect the segment.
 588  */
 589 static int
 590 seg_in_region (const BoxType * b, void *cl)
 591 {
 592   struct info *i = (struct info *) cl;
 593   double y1, y2;
 594   /* for zero slope the search is aligned on the axis so it is already pruned */
 595   if (i->m == 0.)
 596     return 1;
 597   y1 = i->m * b->X1 + i->b;
 598   y2 = i->m * b->X2 + i->b;
 599   if (min (y1, y2) >= b->Y2)
 600     return 0;
 601   if (max (y1, y2) < b->Y1)
 602     return 0;
 603   return 1;                     /* might intersect */
 604 }
 605
 606 /*
 607  * seg_in_seg()
 608  * (C) 2006 harry eaton
 609  * This routine checks if the segment in the tree intersect the search segment.
 610  * If it does, the plines are marked as intersected and the point is marked for
 611  * the cvclist. If the point is not already a vertex, a new vertex is inserted
 612  * and the search for intersections starts over at the beginning.
 613  * That is potentially a significant time penalty, but it does solve the snap rounding
 614  * problem. There are efficient algorithms for finding intersections with snap
 615  * rounding, but I don't have time to implement them right now.
 616  */
 617 static int
 618 seg_in_seg (const BoxType * b, void *cl)
 619 {
 620   struct info *i = (struct info *) cl;
 621   struct seg *s = (struct seg *) b;
 622   Vector s1, s2;
 623   int cnt, res;
 624
 625   cnt = vect_inters2 (s->v->point, s->v->next->point,
 626                       i->v->point, i->v->next->point, s1, s2);
 627   if (!cnt)
 628     return 0;
 629   if (i->touch)                 /* if checking touches one find and we're done */
 630     longjmp (*i->touch, TOUCHES);
 631   i->s->p->Flags.status = ISECTED;
 632   s->p->Flags.status = ISECTED;
 633   for (; cnt; cnt--)
 634     {
 635       res = node_add_point (i->v, s->v, cnt > 1 ? s2 : s1);
 636       if (res < 0)
 637         return 1;               /* error */
 638       /* adjust the bounding box and tree if necessary */
 639       if (res & 2)
 640         {
 641           cntrbox_adjust (i->s->p, cnt > 1 ? s2 : s1);
 642           if (adjust_tree (i->s->p->tree, i->s))
 643             return 1;
 644         }
 645       /* if we added a node in the tree we need to change the tree */
 646       if (res & 1)
 647         {
 648           cntrbox_adjust (s->p, cnt > 1 ? s2 : s1);
 649           if (adjust_tree (i->tree, s))
 650             return 1;
 651         }
 652       if (res & 3)              /* if a point was inserted start over */
 653         {
 654 #ifdef DEBUG_INTERSECT
 655           DEBUGP ("new intersection at (%d, %d)\n", cnt > 1 ? s2[0] : s1[0],
 656                   cnt > 1 ? s2[1] : s1[1]);
 657 #endif
 658           longjmp (i->env, 1);
 659         }
 660     }
 661   return 0;
 662 }
 663
 664 static void *
 665 make_edge_tree (PLINE * pb)
 666 {
 667   struct seg *s;
 668   VNODE *bv;
 669   rtree_t *ans = r_create_tree (NULL, 0, 0);
 670   bv = &pb->head;
 671   do
 672     {
 673       s = malloc (sizeof (struct seg));
 674       if (bv->point[0] < bv->next->point[0])
 675         {
 676           s->box.X1 = bv->point[0];
 677           s->box.X2 = bv->next->point[0] + 1;
 678         }
 679       else
 680         {
 681           s->box.X2 = bv->point[0] + 1;
 682           s->box.X1 = bv->next->point[0];
 683         }
 684       if (bv->point[1] < bv->next->point[1])
 685         {
 686           s->box.Y1 = bv->point[1];
 687           s->box.Y2 = bv->next->point[1] + 1;
 688         }
 689       else
 690         {
 691           s->box.Y2 = bv->point[1] + 1;
 692           s->box.Y1 = bv->next->point[1];
 693         }
 694       s->v = bv;
 695       s->p = pb;
 696       r_insert_entry (ans, (const BoxType *) s, 1);
 697     }
 698   while ((bv = bv->next) != &pb->head);
 699   return (void *) ans;
 700 }
 701
 702 static int
 703 get_seg (const BoxType * b, void *cl)
 704 {
 705   struct info *i = (struct info *) cl;
 706   struct seg *s = (struct seg *) b;
 707   if (i->v == s->v)
 708     {
 709       i->s = s;
 710       longjmp (i->sego, 1);
 711     }
 712   return 0;
 713 }
 714
 715 /*
 716  * intersect()
 717  * (C) 2006, harry eaton
 718  * This uses an rtree to find A-B intersections. Whenever a new vertex is
 719  * added, the search for intersections is re-started because the rounding
 720  * could alter the topology otherwise.
 721  * This should use a faster algorithm for snap rounding intersection finding.
 722  * The best algorthim is probably found in:
 723  *
 724  * "Improved output-sensitive snap rounding," John Hershberger, Proceedings
 725  * of the 22nd annual symposium on Computational geomerty, 2006, pp 357-366.
 726  * http://doi.acm.org/10.1145/1137856.1137909
 727  *
 728  * Algorithms described by de Berg, or Goodrich or Halperin, or Hobby would
 729  * probably work as well.
 730  *
 731  */
 732
 733 static int
 734 intersect (jmp_buf * jb, POLYAREA * b, POLYAREA * a, int add)
 735 {
 736   PLINE *pa, *pb;               /* pline iterators */
 737   PLINE *rtree_over;
 738   PLINE *looping_over;
 739   VNODE *av;                    /* node iterators */
 740   struct info info;
 741   BoxType box;
 742
 743   if (add)
 744     info.touch = NULL;
 745   else
 746     info.touch = jb;
 747   setjmp (info.env);            /* we loop back here whenever a vertex is inserted */
 748   {
 749     pa = a->contours;
 750     pb = b->contours;
 751     while (pa)     /* Loop over the contours of POLYAREA "a" */
 752       {
 753         int found_overlapping_a_b_contour = FALSE;
 754
 755         while (pb) /* Loop over the contours of POLYAREA "b" */
 756           {
 757             /* Are there overlapping bounds? */
 758             if (pb->xmin <= pa->xmax && pb->xmax >= pa->xmin &&
 759                 pb->ymin <= pa->ymax && pb->ymax >= pa->ymin)
 760               {
 761                 found_overlapping_a_b_contour = TRUE;
 762                 break;
 763               }
 764             pb = pb->next;
 765           }
 766
 767         /* If we didn't find anything intersting, move onto the next "a" contour */
 768         if (!found_overlapping_a_b_contour)
 769           {
 770             pa = pa->next;
 771             pb = b->contours;
 772             continue;
 773           }
 774
 775         /* something intersects so check the edges of the contour */
 776
 777         /* Pick which contour has the fewer points, and do the loop
 778          * over that. The r_tree makes hit-testing against a contour
 779          * faster, so we want to do that on the bigger contour.
 780          */
 781         if (pa->Count < pb->Count)
 782           {
 783             rtree_over   = pb;
 784             looping_over = pa;
 785           }
 786         else
 787           {
 788             rtree_over   = pa;
 789             looping_over = pb;
 790           }
 791
 792         av = &looping_over->head;
 793         do  /* Loop over the nodes in the smaller contour */
 794           {
 795             /* check this edge for any insertions */
 796             double dx;
 797             info.v = av;
 798             /* compute the slant for region trimming */
 799             dx = av->next->point[0] - av->point[0];
 800             if (dx == 0)
 801               info.m = 0;
 802             else
 803               {
 804                 info.m = (av->next->point[1] - av->point[1]) / dx;
 805                 info.b = av->point[1] - info.m * av->point[0];
 806               }
 807             box.X2 = (box.X1 = av->point[0]) + 1;
 808             box.Y2 = (box.Y1 = av->point[1]) + 1;
 809
 810             /* fill in the segment in info corresponding to this node */
 811             if (setjmp (info.sego) == 0)
 812               {
 813                 r_search (looping_over->tree, &box, NULL, get_seg, &info);
 814                 assert (0);
 815               }
 816
 817               /* NB: If this actually hits anything, we are teleported back to the beginning */
 818               info.tree = rtree_over->tree;
 819               if (info.tree)
 820                 if (UNLIKELY (r_search (info.tree, &info.s->box,
 821                                         seg_in_region, seg_in_seg, &info)))
 822                   return err_no_memory; /* error */
 823           }
 824         while ((av = av->next) != &looping_over->head);
 825
 826         /* Continue the with the _same_ "a" contour,
 827          * testing it against the next "b" contour.
 828          */
 829         pb = pb->next;
 830       }
 831   }                             /* end of setjmp loop */
 832   return 0;
 833 }
 834
 835 static void
 836 M_POLYAREA_intersect (jmp_buf * e, POLYAREA * afst, POLYAREA * bfst, int add)
 837 {
 838   POLYAREA *a = afst, *b = bfst;
 839   PLINE *curcA, *curcB;
 840   CVCList *the_list = NULL;
 841
 842   if (a == NULL || b == NULL)
 843     error (err_bad_parm);
 844   do
 845     {
 846       do
 847         {
 848           if (a->contours->xmax >= b->contours->xmin &&
 849               a->contours->ymax >= b->contours->ymin &&
 850               a->contours->xmin <= b->contours->xmax &&
 851               a->contours->ymin <= b->contours->ymax)
 852             {
 853               if (UNLIKELY (intersect (e, a, b, add)))
 854                 error (err_no_memory);
 855             }
 856         }
 857       while (add && (a = a->f) != afst);
 858       for (curcB = b->contours; curcB != NULL; curcB = curcB->next)
 859         if (curcB->Flags.status == ISECTED)
 860           {
 861             the_list = add_descriptors (curcB, 'B', the_list);
 862             if (UNLIKELY (the_list == NULL))
 863               error (err_no_memory);
 864           }
 865     }
 866   while (add && (b = b->f) != bfst);
 867   do
 868     {
 869       for (curcA = a->contours; curcA != NULL; curcA = curcA->next)
 870         if (curcA->Flags.status == ISECTED)
 871           {
 872             the_list = add_descriptors (curcA, 'A', the_list);
 873             if (UNLIKELY (the_list == NULL))
 874               error (err_no_memory);
 875           }
 876     }
 877   while (add && (a = a->f) != afst);
 878 }                               /* M_POLYAREA_intersect */
 879
 880 static inline int
 881 cntrbox_inside (PLINE * c1, PLINE * c2)
 882 {
 883   assert (c1 != NULL && c2 != NULL);
 884   return ((c1->xmin >= c2->xmin) &&
 885           (c1->ymin >= c2->ymin) &&
 886           (c1->xmax <= c2->xmax) && (c1->ymax <= c2->ymax));
 887 }
 888
 889 /*****************************************************************/
 890 /* Routines for making labels */
 891
 892 /* cntr_in_M_POLYAREA
 893 returns poly is inside outfst ? TRUE : FALSE */
 894 static int
 895 cntr_in_M_POLYAREA (PLINE * poly, POLYAREA * outfst, BOOLp test)
 896 {
 897   PLINE *curc;
 898   POLYAREA *outer = outfst;
 899   heap_t *heap;
 900
 901   assert (poly != NULL);
 902   assert (outer != NULL);
 903
 904   heap = heap_create ();
 905   do
 906     {
 907       if (cntrbox_inside (poly, outer->contours))
 908         heap_insert (heap, outer->contours->area, (void *) outer);
 909     }
 910   /* if checking touching, use only the first polygon */
 911   while (!test && (outer = outer->f) != outfst);
 912   /* we need only check the smallest poly container
 913    * but we must loop in case the box containter is not
 914    * the poly container */
 915   do
 916     {
 917       if (heap_is_empty (heap))
 918         break;
 919       outer = (POLYAREA *) heap_remove_smallest (heap);
 920       if (poly_ContourInContour (outer->contours, poly))
 921         {
 922           for (curc = outer->contours->next; curc != NULL; curc = curc->next)
 923             if (poly_ContourInContour (curc, poly))
 924               {
 925                 /* it's inside a hole in the smallest polygon
 926                  * no need to check the other polygons */
 927                 heap_destroy (&heap);
 928                 return FALSE;
 929               }
 930           heap_destroy (&heap);
 931           return TRUE;
 932         }
 933     }
 934   while (1);
 935   heap_destroy (&heap);
 936   return FALSE;
 937 }                               /* cntr_in_M_POLYAREA */
 938
 939 #ifdef DEBUG
 940
 941 static char *
 942 theState (VNODE * v)
 943 {
 944   static char u[] = "UNKNOWN";
 945   static char i[] = "INSIDE";
 946   static char o[] = "OUTSIDE";
 947   static char s[] = "SHARED";
 948   static char s2[] = "SHARED2";
 949
 950   switch (NODE_LABEL (v))
 951     {
 952     case INSIDE:
 953       return i;
 954     case OUTSIDE:
 955       return o;
 956     case SHARED:
 957       return s;
 958     case SHARED2:
 959       return s2;
 960     default:
 961       return u;
 962     }
 963 }
 964
 965 static void
 966 print_labels (PLINE * a)
 967 {
 968   VNODE *c = &a->head;
 969
 970   do
 971     {
 972       DEBUGP ("(%ld,%ld)->(%ld,%ld) labeled %s\n", c->point[0], c->point[1],
 973               c->next->point[0], c->next->point[1], theState (c));
 974     }
 975   while ((c = c->next) != &a->head);
 976 }
 977 #endif
 978
 979 /*
 980 label_contour
 981  (C) 2006 harry eaton
 982  (C) 1993 Klamer Schutte
 983  (C) 1997 Alexey Nikitin, Michael Leonov
 984 */
 985
 986 static BOOLp
 987 label_contour (PLINE * a)
 988 {
 989   VNODE *cur = &a->head;
 990   int did_label = FALSE, label = UNKNWN;
 991
 992   do
 993     {
 994       if (cur == &a->head)
 995         did_label = FALSE;
 996       if (NODE_LABEL (cur) != UNKNWN)
 997         {
 998           label = NODE_LABEL (cur);
 999           continue;
1000         }
1001       if (cur->cvc_next)        /* examine cross vertex */
1002         {
1003           label = node_label (cur);
1004           did_label = TRUE;
1005         }
1006       else if (label == INSIDE || label == OUTSIDE)
1007         {
1008           LABEL_NODE (cur, label);
1009           did_label = TRUE;
1010         }
1011     }
1012   while ((cur = cur->next) != &a->head || did_label);
1013 #ifdef DEBUG_ALL_LABELS
1014   print_labels (a);
1015   DEBUGP ("\n\n");
1016 #endif
1017   return FALSE;
1018 }                               /* label_contour */
1019
1020 static BOOLp
1021 cntr_label_POLYAREA (PLINE * poly, POLYAREA * ppl, BOOLp test)
1022 {
1023   assert (ppl != NULL && ppl->contours != NULL);
1024   if (poly->Flags.status == ISECTED)
1025     {
1026       label_contour (poly);     /* should never get here when BOOLp is true */
1027     }
1028   else if (cntr_in_M_POLYAREA (poly, ppl, test))
1029     {
1030       if (test)
1031         return TRUE;
1032       poly->Flags.status = INSIDE;
1033     }
1034   else
1035     {
1036       if (test)
1037         return False;
1038       poly->Flags.status = OUTSIDE;
1039     }
1040   return FALSE;
1041 }                               /* cntr_label_POLYAREA */
1042
1043 static BOOLp
1044 M_POLYAREA_label (POLYAREA * afst, POLYAREA * b, BOOLp touch)
1045 {
1046   POLYAREA *a = afst;
1047   PLINE *curc;
1048
1049   assert (a != NULL);
1050   do
1051     {
1052       for (curc = a->contours; curc != NULL; curc = curc->next)
1053         if (cntr_label_POLYAREA (curc, b, touch))
1054           {
1055             if (touch)
1056               return TRUE;
1057           }
1058     }
1059   while (!touch && (a = a->f) != afst);
1060   return FALSE;
1061 }
1062
1063 /****************************************************************/
1064
1065 /* routines for temporary storing resulting contours */
1066 static void
1067 InsCntr (jmp_buf * e, PLINE * c, POLYAREA ** dst)
1068 {
1069   POLYAREA *newp;
1070
1071   if (*dst == NULL)
1072     {
1073       MemGet (*dst, POLYAREA);
1074       (*dst)->f = (*dst)->b = *dst;
1075       newp = *dst;
1076     }
1077   else
1078     {
1079       MemGet (newp, POLYAREA);
1080       newp->f = *dst;
1081       newp->b = (*dst)->b;
1082       newp->f->b = newp->b->f = newp;
1083     }
1084   newp->contours = c;
1085   c->next = NULL;
1086 }                               /* InsCntr */
1087
1088 static void
1089 PutContour (jmp_buf * e, PLINE * cntr, POLYAREA ** contours, PLINE ** holes,
1090             PLINE * parent)
1091 {
1092   assert (cntr != NULL);
1093   assert (cntr->Count > 2);
1094   cntr->next = NULL;
1095   if (cntr->Flags.orient == PLF_DIR)
1096     InsCntr (e, cntr, contours);
1097   /* put hole into temporary list */
1098   else
1099     {
1100       /* if we know this belongs inside the parent, put it there now */
1101       if (parent)
1102         {
1103           cntr->next = parent->next;
1104           parent->next = cntr;
1105         }
1106       else
1107         {
1108           cntr->next = *holes;
1109           *holes = cntr;        /* let cntr be 1st hole in list */
1110         }
1111     }
1112 }                               /* PutContour */
1113
1114 static int
1115 heap_it (const BoxType * b, void *cl)
1116 {
1117   heap_t *heap = (heap_t *) cl;
1118   PLINE *p = (PLINE *) b;
1119   if (p->Count == 0)
1120     return 0;  /* how did this happen? */
1121   heap_insert (heap, p->area, (void *) p);
1122   return 1;
1123 }
1124
1125 static void
1126 InsertHoles (jmp_buf * e, POLYAREA * dest, PLINE ** src)
1127 {
1128   POLYAREA *curc;
1129   PLINE *curh, *container, *tmp;
1130   heap_t *heap;
1131   rtree_t *tree;
1132
1133   if (*src == NULL)
1134     return;                     /* empty hole list */
1135   if (dest == NULL)
1136     error (err_bad_parm);       /* empty contour list */
1137
1138   /* make an rtree of contours */
1139   tree = r_create_tree (NULL, 0, 0);
1140   curc = dest;
1141   do
1142     {
1143       r_insert_entry (tree, (const BoxType *) curc->contours, 0);
1144     }
1145   while ((curc = curc->f) != dest);
1146   /* loop through the holes and put them where they belong */
1147   while ((curh = *src) != NULL)
1148     {
1149       *src = curh->next;
1150
1151       container = NULL;
1152       /* build a heap of all of the polys that the hole is inside its bounding box */
1153       heap = heap_create ();
1154       r_search (tree, (BoxType *) curh, NULL, heap_it, heap);
1155       if (heap_is_empty (heap))
1156         {
1157 #ifndef NDEBUG
1158 #ifdef DEBUG
1159           poly_dump (dest);
1160 #endif
1161 #endif
1162           poly_DelContour (&curh);
1163           error (err_bad_parm);
1164         }
1165       /* Now search the heap for the container. If there was only one item
1166        * in the heap, assume it is the container without the expense of
1167        * proving it.
1168        */
1169       tmp = (PLINE *) heap_remove_smallest (heap);
1170       if (heap_is_empty (heap))
1171         {                       /* only one possibility it must be the right one */
1172           assert (poly_ContourInContour (tmp, curh));
1173           container = tmp;
1174         }
1175       else
1176         {
1177           do
1178             {
1179               if (poly_ContourInContour (tmp, curh))
1180                 {
1181                   container = tmp;
1182                   break;
1183                 }
1184               if (heap_is_empty (heap))
1185                 break;
1186               tmp = (PLINE *) heap_remove_smallest (heap);
1187             }
1188           while (1);
1189         }
1190       heap_destroy (&heap);
1191       if (container == NULL)
1192         {
1193           /* bad input polygons were given */
1194 #ifndef NDEBUG
1195 #ifdef DEBUG
1196           poly_dump (dest);
1197 #endif
1198 #endif
1199           curh->next = NULL;
1200           poly_DelContour (&curh);
1201           error (err_bad_parm);
1202         }
1203       else
1204         {
1205           /* link at front of hole list */
1206           tmp = container->next;
1207           container->next = curh;
1208           curh->next = tmp;
1209         }
1210     }
1211   r_destroy_tree (&tree);
1212 }                               /* InsertHoles */
1213
1214
1215 /****************************************************************/
1216 /* routines for collecting result */
1217
1218 typedef enum
1219 {
1220   FORW, BACKW
1221 } DIRECTION;
1222
1223 /* Start Rule */
1224 typedef int (*S_Rule) (VNODE *, DIRECTION *);
1225
1226 /* Jump Rule  */
1227 typedef int (*J_Rule) (char, VNODE *, DIRECTION *);
1228
1229 static int
1230 UniteS_Rule (VNODE * cur, DIRECTION * initdir)
1231 {
1232   *initdir = FORW;
1233   return (NODE_LABEL (cur) == OUTSIDE) || (NODE_LABEL (cur) == SHARED);
1234 }
1235
1236 static int
1237 IsectS_Rule (VNODE * cur, DIRECTION * initdir)
1238 {
1239   *initdir = FORW;
1240   return (NODE_LABEL (cur) == INSIDE) || (NODE_LABEL (cur) == SHARED);
1241 }
1242
1243 static int
1244 SubS_Rule (VNODE * cur, DIRECTION * initdir)
1245 {
1246   *initdir = FORW;
1247   return (NODE_LABEL (cur) == OUTSIDE) || (NODE_LABEL (cur) == SHARED2);
1248 }
1249
1250 static int
1251 XorS_Rule (VNODE * cur, DIRECTION * initdir)
1252 {
1253   if (cur->Flags.status == INSIDE)
1254     {
1255       *initdir = BACKW;
1256       return TRUE;
1257     }
1258   if (cur->Flags.status == OUTSIDE)
1259     {
1260       *initdir = FORW;
1261       return TRUE;
1262     }
1263   return FALSE;
1264 }
1265
1266 static int
1267 IsectJ_Rule (char p, VNODE * v, DIRECTION * cdir)
1268 {
1269   assert (*cdir == FORW);
1270   return (v->Flags.status == INSIDE || v->Flags.status == SHARED);
1271 }
1272
1273 static int
1274 UniteJ_Rule (char p, VNODE * v, DIRECTION * cdir)
1275 {
1276   assert (*cdir == FORW);
1277   return (v->Flags.status == OUTSIDE || v->Flags.status == SHARED);
1278 }
1279
1280 static int
1281 XorJ_Rule (char p, VNODE * v, DIRECTION * cdir)
1282 {
1283   if (v->Flags.status == INSIDE)
1284     {
1285       *cdir = BACKW;
1286       return TRUE;
1287     }
1288   if (v->Flags.status == OUTSIDE)
1289     {
1290       *cdir = FORW;
1291       return TRUE;
1292     }
1293   return FALSE;
1294 }
1295
1296 static int
1297 SubJ_Rule (char p, VNODE * v, DIRECTION * cdir)
1298 {
1299   if (p == 'B' && v->Flags.status == INSIDE)
1300     {
1301       *cdir = BACKW;
1302       return TRUE;
1303     }
1304   if (p == 'A' && v->Flags.status == OUTSIDE)
1305     {
1306       *cdir = FORW;
1307       return TRUE;
1308     }
1309   if (v->Flags.status == SHARED2)
1310     {
1311       if (p == 'A')
1312         *cdir = FORW;
1313       else
1314         *cdir = BACKW;
1315       return TRUE;
1316     }
1317   return FALSE;
1318 }
1319
1320 /* return the edge that comes next.
1321  * if the direction is BACKW, then we return the next vertex
1322  * so that prev vertex has the flags for the edge
1323  *
1324  * returns true if an edge is found, false otherwise
1325  */
1326 static int
1327 jump (VNODE ** cur, DIRECTION * cdir, J_Rule rule)
1328 {
1329   CVCList *d, *start;
1330   VNODE *e;
1331   DIRECTION new;
1332
1333   if (!(*cur)->cvc_prev)        /* not a cross-vertex */
1334     {
1335       if (*cdir == FORW ? (*cur)->Flags.mark : (*cur)->prev->Flags.mark)
1336         return FALSE;
1337       return TRUE;
1338     }
1339 #ifdef DEBUG_JUMP
1340   DEBUGP ("jump entering node at (%ld, %ld)\n", (*cur)->point[0],
1341           (*cur)->point[1]);
1342 #endif
1343   if (*cdir == FORW)
1344     d = (*cur)->cvc_prev->prev;
1345   else
1346     d = (*cur)->cvc_next->prev;
1347   start = d;
1348   do
1349     {
1350       e = d->parent;
1351       if (d->side == 'P')
1352         e = e->prev;
1353       new = *cdir;
1354       if (!e->Flags.mark && rule (d->poly, e, &new))
1355         {
1356           if ((d->side == 'N' && new == FORW) ||
1357               (d->side == 'P' && new == BACKW))
1358             {
1359 #ifdef DEBUG_JUMP
1360               if (new == FORW)
1361                 DEBUGP ("jump leaving node at (%ld, %ld)\n",
1362                         e->next->point[0], e->next->point[1]);
1363               else
1364                 DEBUGP ("jump leaving node at (%ld, %ld)\n",
1365                         e->point[0], e->point[1]);
1366 #endif
1367               *cur = d->parent;
1368               *cdir = new;
1369               return TRUE;
1370             }
1371         }
1372     }
1373   while ((d = d->prev) != start);
1374   return FALSE;
1375 }
1376
1377 static int
1378 Gather (VNODE * start, PLINE ** result, J_Rule v_rule, DIRECTION initdir)
1379 {
1380   VNODE *cur = start, *newn;
1381   DIRECTION dir = initdir;
1382 #ifdef DEBUG_GATHER
1383   DEBUGP ("gather direction = %d\n", dir);
1384 #endif
1385   assert (*result == NULL);
1386   do
1387     {
1388       /* see where to go next */
1389       if (!jump (&cur, &dir, v_rule))
1390         break;
1391       /* add edge to polygon */
1392       if (!*result)
1393         {
1394           *result = poly_NewContour (cur->point);
1395           if (*result == NULL)
1396             return err_no_memory;
1397         }
1398       else
1399         {
1400           if ((newn = poly_CreateNode (cur->point)) == NULL)
1401             return err_no_memory;
1402           poly_InclVertex ((*result)->head.prev, newn);
1403         }
1404 #ifdef DEBUG_GATHER
1405       DEBUGP ("gather vertex at (%ld, %ld)\n", cur->point[0], cur->point[1]);
1406 #endif
1407       /* Now mark the edge as included.  */
1408       newn = (dir == FORW ? cur : cur->prev);
1409       newn->Flags.mark = 1;
1410       /* for SHARED edge mark both */
1411       if (newn->shared)
1412         newn->shared->Flags.mark = 1;
1413
1414       /* Advance to the next edge.  */
1415       cur = (dir == FORW ? cur->next : cur->prev);
1416     }
1417   while (1);
1418   return err_ok;
1419 }                               /* Gather */
1420
1421 static void
1422 Collect1 (jmp_buf * e, VNODE *cur, DIRECTION dir, POLYAREA **contours, PLINE ** holes, J_Rule j_rule)
1423 {
1424   PLINE *p = NULL;              /* start making contour */
1425   int errc = err_ok;
1426         if ((errc =
1427              Gather (dir == FORW ? cur : cur->next, &p, j_rule,
1428                      dir)) != err_ok)
1429           {
1430             if (p != NULL)
1431               poly_DelContour (&p);
1432             error (errc);
1433           }
1434         if (!p)
1435           return;
1436         poly_PreContour (p, TRUE);
1437         if (p->Count > 2)
1438           {
1439 #ifdef DEBUG_GATHER
1440             DEBUGP ("adding contour with %d verticies and direction %c\n",
1441                     p->Count, p->Flags.orient ? 'F' : 'B');
1442 #endif
1443             PutContour (e, p, contours, holes, NULL);
1444           }
1445         else
1446           {
1447             /* some sort of computation error ? */
1448 #ifdef DEBUG_GATHER
1449             DEBUGP ("Bad contour! Not enough points!!\n");
1450 #endif
1451             poly_DelContour (&p);
1452           }
1453 }
1454
1455 static void
1456 Collect (jmp_buf * e, PLINE * a, POLYAREA ** contours, PLINE ** holes,
1457          S_Rule s_rule, J_Rule j_rule)
1458 {
1459   VNODE *cur, *other;
1460   DIRECTION dir;
1461
1462   cur = &a->head;
1463   do
1464    {
1465     if (s_rule (cur, &dir) && cur->Flags.mark == 0)
1466         Collect1(e, cur, dir, contours, holes, j_rule);
1467     other = cur;
1468     if ((other->cvc_prev && jump(&other, &dir, j_rule)))
1469         Collect1(e, other, dir, contours, holes, j_rule);
1470    }
1471   while ((cur = cur->next) != &a->head);
1472 }                               /* Collect */
1473
1474
1475 static int
1476 cntr_Collect (jmp_buf * e, PLINE ** A, POLYAREA ** contours, PLINE ** holes,
1477               int action, PLINE * parent)
1478 {
1479   PLINE *tmprev;
1480
1481   if ((*A)->Flags.status == ISECTED)
1482     {
1483       switch (action)
1484         {
1485         case PBO_UNITE:
1486           Collect (e, *A, contours, holes, UniteS_Rule, UniteJ_Rule);
1487           break;
1488         case PBO_ISECT:
1489           Collect (e, *A, contours, holes, IsectS_Rule, IsectJ_Rule);
1490           break;
1491         case PBO_XOR:
1492           Collect (e, *A, contours, holes, XorS_Rule, XorJ_Rule);
1493           break;
1494         case PBO_SUB:
1495           Collect (e, *A, contours, holes, SubS_Rule, SubJ_Rule);
1496           break;
1497         };
1498     }
1499   else
1500     {
1501       switch (action)
1502         {
1503         case PBO_ISECT:
1504           if ((*A)->Flags.status == INSIDE)
1505             {
1506               tmprev = *A;
1507               /* disappear this contour */
1508               *A = tmprev->next;
1509               tmprev->next = NULL;
1510               PutContour (e, tmprev, contours, holes, NULL);
1511               return TRUE;
1512             }
1513           break;
1514         case PBO_XOR:
1515           if ((*A)->Flags.status == INSIDE)
1516             {
1517               tmprev = *A;
1518               /* disappear this contour */
1519               *A = tmprev->next;
1520               tmprev->next = NULL;
1521               poly_InvContour (tmprev);
1522               PutContour (e, tmprev, contours, holes, NULL);
1523               return TRUE;
1524             }
1525           break;
1526         case PBO_UNITE:
1527         case PBO_SUB:
1528           if ((*A)->Flags.status == OUTSIDE)
1529             {
1530               tmprev = *A;
1531               /* disappear this contour */
1532               *A = tmprev->next;
1533               tmprev->next = NULL;
1534               PutContour (e, tmprev, contours, holes, parent);
1535               return TRUE;
1536             }
1537           break;
1538         }
1539     }
1540   return FALSE;
1541 }                               /* cntr_Collect */
1542
1543 static void
1544 M_B_AREA_Collect (jmp_buf * e, POLYAREA * bfst, POLYAREA ** contours,
1545                   PLINE ** holes, int action)
1546 {
1547   POLYAREA *b = bfst;
1548   PLINE **cur, **next, *tmp;
1549
1550   assert (b != NULL);
1551   do
1552     {
1553       for (cur = &b->contours; *cur != NULL; cur = next)
1554         {
1555           next = &((*cur)->next);
1556           if ((*cur)->Flags.status == ISECTED)
1557             continue;
1558
1559           if ((*cur)->Flags.status == INSIDE)
1560             switch (action)
1561               {
1562               case PBO_XOR:
1563               case PBO_SUB:
1564                 poly_InvContour (*cur);
1565               case PBO_ISECT:
1566                 tmp = *cur;
1567                 *cur = tmp->next;
1568                 next = cur;
1569                 tmp->next = NULL;
1570                 tmp->Flags.status = UNKNWN;
1571                 PutContour (e, tmp, contours, holes, NULL);
1572                 break;
1573               case PBO_UNITE:
1574                 break;          /* nothing to do - already included */
1575               }
1576           else if ((*cur)->Flags.status == OUTSIDE)
1577             switch (action)
1578               {
1579               case PBO_XOR:
1580               case PBO_UNITE:
1581                 /* include */
1582                 tmp = *cur;
1583                 *cur = tmp->next;
1584                 next = cur;
1585                 tmp->next = NULL;
1586                 tmp->Flags.status = UNKNWN;
1587                 PutContour (e, tmp, contours, holes, NULL);
1588                 break;
1589               case PBO_ISECT:
1590               case PBO_SUB:
1591                 break;          /* do nothing, not included */
1592               }
1593         }
1594     }
1595   while ((b = b->f) != bfst);
1596 }
1597
1598
1599 static void
1600 M_POLYAREA_Collect (jmp_buf * e, POLYAREA * afst, POLYAREA ** contours,
1601                     PLINE ** holes, int action, BOOLp maybe)
1602 {
1603   POLYAREA *a = afst;
1604   PLINE **cur, **next, *parent;
1605
1606   assert (a != NULL);
1607   while ((a = a->f) != afst);
1608   /* now the non-intersect parts are collected in temp/holes */
1609   do
1610     {
1611       if (maybe && a->contours->Flags.status != ISECTED)
1612         parent = a->contours;
1613       else
1614         parent = NULL;
1615       for (cur = &a->contours; *cur != NULL; cur = next)
1616         {
1617           next = &((*cur)->next);
1618           /* if we disappear a contour, don't advance twice */
1619           if (cntr_Collect
1620               (e, cur, contours, holes, action,
1621                *cur == parent ? NULL : parent))
1622             next = cur;
1623         }
1624     }
1625   while ((a = a->f) != afst);
1626 }
1627
1628 /* determine if two polygons touch or overlap */
1629 BOOLp
1630 Touching (POLYAREA * a, POLYAREA * b)
1631 {
1632   jmp_buf e;
1633   int code;
1634
1635   if ((code = setjmp (e)) == 0)
1636     {
1637 #ifdef DEBUG
1638       if (!poly_Valid (a))
1639         return -1;
1640       if (!poly_Valid (b))
1641         return -1;
1642 #endif
1643       M_POLYAREA_intersect (&e, a, b, False);
1644
1645       if (M_POLYAREA_label (a, b, TRUE))
1646         return TRUE;
1647       if (M_POLYAREA_label (b, a, TRUE))
1648         return TRUE;
1649     }
1650   else if (code == TOUCHES)
1651     return TRUE;
1652   return FALSE;
1653 }
1654
1655 /* the main clipping routines */
1656 int
1657 poly_Boolean (const POLYAREA * a_org, const POLYAREA * b_org,
1658               POLYAREA ** res, int action)
1659 {
1660   POLYAREA *a = NULL, *b = NULL;
1661
1662   if (!poly_M_Copy0 (&a, a_org) || !poly_M_Copy0 (&b, b_org))
1663     return err_no_memory;
1664
1665   return poly_Boolean_free (a, b, res, action);
1666 }                               /* poly_Boolean */
1667
1668 /* just like poly_Boolean but frees the input polys */
1669 int
1670 poly_Boolean_free (POLYAREA * ai, POLYAREA * bi, POLYAREA ** res, int action)
1671 {
1672   POLYAREA *a = ai, *b = bi;
1673   PLINE *p, *holes = NULL;
1674   jmp_buf e;
1675   int code;
1676
1677   *res = NULL;
1678
1679   if (!a)
1680     {
1681       switch (action)
1682         {
1683         case PBO_XOR:
1684         case PBO_UNITE:
1685           *res = bi;
1686           return err_ok;
1687         case PBO_SUB:
1688         case PBO_ISECT:
1689           if (b != NULL)
1690             poly_Free (&b);
1691           return err_ok;
1692         }
1693     }
1694   if (!b)
1695     {
1696       switch (action)
1697         {
1698         case PBO_SUB:
1699         case PBO_XOR:
1700         case PBO_UNITE:
1701           *res = ai;
1702           return err_ok;
1703         case PBO_ISECT:
1704           if (a != NULL)
1705             poly_Free (&a);
1706           return err_ok;
1707         }
1708     }
1709
1710   if ((code = setjmp (e)) == 0)
1711     {
1712 #ifdef DEBUG
1713       assert (poly_Valid (a));
1714       assert (poly_Valid (b));
1715 #endif
1716
1717       M_POLYAREA_intersect (&e, a, b, TRUE);
1718
1719       M_POLYAREA_label (a, b, FALSE);
1720       M_POLYAREA_label (b, a, FALSE);
1721
1722       M_POLYAREA_Collect (&e, a, res, &holes, action, b->f == b
1723                           && !b->contours->next
1724                           && b->contours->Flags.status != ISECTED);
1725       poly_Free (&a);
1726       M_B_AREA_Collect (&e, b, res, &holes, action);
1727       poly_Free (&b);
1728
1729       InsertHoles (&e, *res, &holes);
1730     }
1731   /* delete holes if any left */
1732   while ((p = holes) != NULL)
1733     {
1734       holes = p->next;
1735       poly_DelContour (&p);
1736     }
1737
1738   if (code)
1739     {
1740       poly_Free (res);
1741       return code;
1742     }
1743   assert (!*res || poly_Valid (*res));
1744   return code;
1745 }                               /* poly_Boolean_free */
1746
1747 static void
1748 clear_marks (POLYAREA * p)
1749 {
1750   POLYAREA *n = p;
1751   PLINE *c;
1752   VNODE *v;
1753
1754   do
1755     {
1756       for (c = n->contours; c; c = c->next)
1757         {
1758           v = &c->head;
1759           do
1760             {
1761               v->Flags.mark = 0;
1762             }
1763           while ((v = v->next) != &c->head);
1764         }
1765     }
1766   while ((n = n->f) != p);
1767 }
1768
1769 /* compute the intersection and subtraction (divides "a" into two pieces)
1770  * and frees the input polys. This assumes that bi is a single simple polygon.
1771  */
1772 int
1773 poly_AndSubtract_free (POLYAREA * ai, POLYAREA * bi,
1774                        POLYAREA ** aandb, POLYAREA ** aminusb)
1775 {
1776   POLYAREA *a = ai, *b = bi;
1777   PLINE *p, *holes = NULL;
1778   jmp_buf e;
1779   int code;
1780
1781   *aandb = NULL;
1782   *aminusb = NULL;
1783
1784   if ((code = setjmp (e)) == 0)
1785     {
1786
1787 #ifdef DEBUG
1788       if (!poly_Valid (a))
1789         return -1;
1790       if (!poly_Valid (b))
1791         return -1;
1792 #endif
1793       M_POLYAREA_intersect (&e, a, b, TRUE);
1794
1795       M_POLYAREA_label (a, b, FALSE);
1796       M_POLYAREA_label (b, a, FALSE);
1797
1798       M_POLYAREA_Collect (&e, a, aandb, &holes, PBO_ISECT, FALSE);
1799       InsertHoles (&e, *aandb, &holes);
1800       assert (poly_Valid (*aandb));
1801       /* delete holes if any left */
1802       while ((p = holes) != NULL)
1803         {
1804           holes = p->next;
1805           poly_DelContour (&p);
1806         }
1807       holes = NULL;
1808       clear_marks (a);
1809       clear_marks (b);
1810       M_POLYAREA_Collect (&e, a, aminusb, &holes, PBO_SUB, FALSE);
1811       InsertHoles (&e, *aminusb, &holes);
1812       poly_Free (&a);
1813       poly_Free (&b);
1814       assert (poly_Valid (*aminusb));
1815     }
1816   /* delete holes if any left */
1817   while ((p = holes) != NULL)
1818     {
1819       holes = p->next;
1820       poly_DelContour (&p);
1821     }
1822
1823
1824   if (code)
1825     {
1826       poly_Free (aandb);
1827       poly_Free (aminusb);
1828       return code;
1829     }
1830   assert (!*aandb || poly_Valid (*aandb));
1831   assert (!*aminusb || poly_Valid (*aminusb));
1832   return code;
1833 }                               /* poly_AndSubtract_free */
1834
1835 static inline int
1836 cntrbox_pointin (PLINE * c, Vector p)
1837 {
1838   return (p[0] >= c->xmin && p[1] >= c->ymin &&
1839           p[0] <= c->xmax && p[1] <= c->ymax);
1840
1841 }
1842
1843 static inline int
1844 node_neighbours (VNODE * a, VNODE * b)
1845 {
1846   return (a == b) || (a->next == b) || (b->next == a) || (a->next == b->next);
1847 }
1848
1849 VNODE *
1850 poly_CreateNode (Vector v)
1851 {
1852   VNODE *res;
1853   register int *c;
1854
1855   assert (v);
1856   res = (VNODE *) calloc (1, sizeof (VNODE));
1857   if (res == NULL)
1858     return NULL;
1859   // bzero (res, sizeof (VNODE) - sizeof(Vector));
1860   c = res->point;
1861   *c++ = *v++;
1862   *c = *v;
1863   return res;
1864 }
1865
1866 void
1867 poly_IniContour (PLINE * c)
1868 {
1869   if (c == NULL)
1870     return;
1871   /* bzero (c, sizeof(PLINE)); */
1872   c->head.next = c->head.prev = &c->head;
1873   c->xmin = c->ymin = 0x7fffffff;
1874   c->xmax = c->ymax = 0x80000000;
1875 }
1876
1877 PLINE *
1878 poly_NewContour (Vector v)
1879 {
1880   PLINE *res;
1881
1882   res = (PLINE *) calloc (1, sizeof (PLINE));
1883   if (res == NULL)
1884     return NULL;
1885
1886   poly_IniContour (res);
1887
1888   if (v != NULL)
1889     {
1890       Vcopy (res->head.point, v);
1891       cntrbox_adjust (res, v);
1892     }
1893
1894   return res;
1895 }
1896
1897 void
1898 poly_ClrContour (PLINE * c)
1899 {
1900   VNODE *cur;
1901
1902   assert (c != NULL);
1903   while ((cur = c->head.next) != &c->head)
1904     {
1905       poly_ExclVertex (cur);
1906       free (cur);
1907     }
1908   poly_IniContour (c);
1909 }
1910
1911 void
1912 poly_DelContour (PLINE ** c)
1913 {
1914   VNODE *cur, *prev;
1915
1916   if (*c == NULL)
1917     return;
1918   for (cur = (*c)->head.prev; cur != &(*c)->head; cur = prev)
1919     {
1920       prev = cur->prev;
1921       if (cur->cvc_next != NULL)
1922         {
1923           free (cur->cvc_next);
1924           free (cur->cvc_prev);
1925         }
1926       free (cur);
1927     }
1928   if ((*c)->head.cvc_next != NULL)
1929     {
1930       free ((*c)->head.cvc_next);
1931       free ((*c)->head.cvc_prev);
1932     }
1933   /* FIXME -- strict aliasing violation.  */
1934   if ((*c)->tree)
1935     {
1936       rtree_t *r = (*c)->tree;
1937       r_destroy_tree (&r);
1938     }
1939   free (*c), *c = NULL;
1940 }
1941
1942 void
1943 poly_PreContour (PLINE * C, BOOLp optimize)
1944 {
1945   double area = 0;
1946   VNODE *p, *c;
1947   Vector p1, p2;
1948
1949   assert (C != NULL);
1950
1951   if (optimize)
1952     {
1953       for (c = (p = &C->head)->next; c != &C->head; c = (p = c)->next)
1954         {
1955           /* if the previous node is on the same line with this one, we should remove it */
1956           Vsub2 (p1, c->point, p->point);
1957           Vsub2 (p2, c->next->point, c->point);
1958           /* If the product below is zero then
1959            * the points on either side of c
1960            * are on the same line!
1961            * So, remove the point c
1962            */
1963
1964           if (vect_det2 (p1, p2) == 0)
1965             {
1966               poly_ExclVertex (c);
1967               free (c);
1968               c = p;
1969             }
1970         }
1971     }
1972   C->Count = 0;
1973   C->xmin = C->xmax = C->head.point[0];
1974   C->ymin = C->ymax = C->head.point[1];
1975
1976   p = (c = &C->head)->prev;
1977   if (c != p)
1978     {
1979       do
1980         {
1981           /* calculate area for orientation */
1982           area +=
1983             (double) (p->point[0] - c->point[0]) * (p->point[1] +
1984                                                     c->point[1]);
1985           cntrbox_adjust (C, c->point);
1986           C->Count++;
1987         }
1988       while ((c = (p = c)->next) != &C->head);
1989     }
1990   C->area = ABS (area);
1991   if (C->Count > 2)
1992     C->Flags.orient = ((area < 0) ? PLF_INV : PLF_DIR);
1993   C->tree = make_edge_tree (C);
1994 }                               /* poly_PreContour */
1995
1996 static int
1997 flip_cb (const BoxType * b, void *cl)
1998 {
1999   struct seg *s = (struct seg *) b;
2000   s->v = s->v->prev;
2001   return 1;
2002 }
2003
2004 void
2005 poly_InvContour (PLINE * c)
2006 {
2007   VNODE *cur, *next;
2008   int r;
2009
2010   assert (c != NULL);
2011   cur = &c->head;
2012   do
2013     {
2014       next = cur->next;
2015       cur->next = cur->prev;
2016       cur->prev = next;
2017       /* fix the segment tree */
2018     }
2019   while ((cur = next) != &c->head);
2020   c->Flags.orient ^= 1;
2021   if (c->tree)
2022     {
2023       r = r_search (c->tree, NULL, NULL, flip_cb, NULL);
2024       assert (r == c->Count);
2025     }
2026 }
2027
2028 void
2029 poly_ExclVertex (VNODE * node)
2030 {
2031   assert (node != NULL);
2032   if (node->cvc_next)
2033     {
2034       free (node->cvc_next);
2035       free (node->cvc_prev);
2036     }
2037   node->prev->next = node->next;
2038   node->next->prev = node->prev;
2039 }
2040
2041 void
2042 poly_InclVertex (VNODE * after, VNODE * node)
2043 {
2044   double a, b;
2045   assert (after != NULL);
2046   assert (node != NULL);
2047
2048   node->prev = after;
2049   node->next = after->next;
2050   after->next = after->next->prev = node;
2051   /* remove points on same line */
2052   if (node->prev->prev == node)
2053     return;                     /* we don't have 3 points in the poly yet */
2054   a = (node->point[1] - node->prev->prev->point[1]);
2055   a *= (node->prev->point[0] - node->prev->prev->point[0]);
2056   b = (node->point[0] - node->prev->prev->point[0]);
2057   b *= (node->prev->point[1] - node->prev->prev->point[1]);
2058   if (fabs (a - b) < EPSILON)
2059     {
2060       VNODE *t = node->prev;
2061       t->prev->next = node;
2062       node->prev = t->prev;
2063       free (t);
2064     }
2065 }
2066
2067 BOOLp
2068 poly_CopyContour (PLINE ** dst, PLINE * src)
2069 {
2070   VNODE *cur, *newnode;
2071
2072   assert (src != NULL);
2073   *dst = NULL;
2074   *dst = poly_NewContour (src->head.point);
2075   if (*dst == NULL)
2076     return FALSE;
2077
2078   (*dst)->Count = src->Count;
2079   (*dst)->Flags.orient = src->Flags.orient;
2080   (*dst)->xmin = src->xmin, (*dst)->xmax = src->xmax;
2081   (*dst)->ymin = src->ymin, (*dst)->ymax = src->ymax;
2082   (*dst)->area = src->area;
2083
2084   for (cur = src->head.next; cur != &src->head; cur = cur->next)
2085     {
2086       if ((newnode = poly_CreateNode (cur->point)) == NULL)
2087         return FALSE;
2088       // newnode->Flags = cur->Flags;
2089       poly_InclVertex ((*dst)->head.prev, newnode);
2090     }
2091   (*dst)->tree = make_edge_tree (*dst);
2092   return TRUE;
2093 }
2094
2095 /**********************************************************************/
2096 /* polygon routines */
2097
2098 BOOLp
2099 poly_Copy0 (POLYAREA ** dst, const POLYAREA * src)
2100 {
2101   *dst = NULL;
2102   if (src != NULL)
2103     *dst = calloc (1, sizeof (POLYAREA));
2104   if (*dst == NULL)
2105     return FALSE;
2106
2107   return poly_Copy1 (*dst, src);
2108 }
2109
2110 BOOLp
2111 poly_Copy1 (POLYAREA * dst, const POLYAREA * src)
2112 {
2113   PLINE *cur, **last = &dst->contours;
2114
2115   *last = NULL;
2116   dst->f = dst->b = dst;
2117
2118   for (cur = src->contours; cur != NULL; cur = cur->next)
2119     {
2120       if (!poly_CopyContour (last, cur))
2121         return FALSE;
2122       last = &(*last)->next;
2123     }
2124   return TRUE;
2125 }
2126
2127 void
2128 poly_M_Incl (POLYAREA ** list, POLYAREA * a)
2129 {
2130   if (*list == NULL)
2131     a->f = a->b = a, *list = a;
2132   else
2133     {
2134       a->f = *list;
2135       a->b = (*list)->b;
2136       (*list)->b = (*list)->b->f = a;
2137     }
2138 }
2139
2140 BOOLp
2141 poly_M_Copy0 (POLYAREA ** dst, const POLYAREA * srcfst)
2142 {
2143   const POLYAREA *src = srcfst;
2144   POLYAREA *di;
2145
2146   *dst = NULL;
2147   if (src == NULL)
2148     return FALSE;
2149   do
2150     {
2151       if ((di = poly_Create ()) == NULL || !poly_Copy1 (di, src))
2152         return FALSE;
2153       poly_M_Incl (dst, di);
2154     }
2155   while ((src = src->f) != srcfst);
2156   return TRUE;
2157 }
2158
2159 BOOLp
2160 poly_InclContour (POLYAREA * p, PLINE * c)
2161 {
2162   PLINE *tmp;
2163
2164   if ((c == NULL) || (p == NULL))
2165     return FALSE;
2166   if (c->Flags.orient == PLF_DIR)
2167     {
2168       if (p->contours != NULL)
2169         return FALSE;
2170       p->contours = c;
2171     }
2172   else
2173     {
2174       if (p->contours == NULL)
2175         return FALSE;
2176       /* link at front of hole list */
2177       tmp = p->contours->next;
2178       p->contours->next = c;
2179       c->next = tmp;
2180     }
2181   return TRUE;
2182 }
2183
2184 typedef struct pip
2185 {
2186   int f;
2187   Vector p;
2188   jmp_buf env;
2189 } pip;
2190
2191
2192 static int
2193 crossing (const BoxType * b, void *cl)
2194 {
2195   struct seg *s = (struct seg *) b;
2196   struct pip *p = (struct pip *) cl;
2197
2198   if (s->v->point[1] <= p->p[1])
2199     {
2200       if (s->v->next->point[1] > p->p[1])
2201         {
2202           Vector v1, v2;
2203           long long cross;
2204           Vsub2 (v1, s->v->next->point, s->v->point);
2205           Vsub2 (v2, p->p, s->v->point);
2206           cross = (long long) v1[0] * v2[1] - (long long) v2[0] * v1[1];
2207           if (cross == 0)
2208             {
2209               p->f = 1;
2210               longjmp (p->env, 1);
2211             }
2212           if (cross > 0)
2213             p->f += 1;
2214         }
2215     }
2216   else
2217     {
2218       if (s->v->next->point[1] <= p->p[1])
2219         {
2220           Vector v1, v2;
2221           long long cross;
2222           Vsub2 (v1, s->v->next->point, s->v->point);
2223           Vsub2 (v2, p->p, s->v->point);
2224           cross = (long long) v1[0] * v2[1] - (long long) v2[0] * v1[1];
2225           if (cross == 0)
2226             {
2227               p->f = 1;
2228               longjmp (p->env, 1);
2229             }
2230           if (cross < 0)
2231             p->f -= 1;
2232         }
2233     }
2234   return 1;
2235 }
2236
2237 int
2238 poly_InsideContour (PLINE * c, Vector p)
2239 {
2240   struct pip info;
2241   BoxType ray;
2242
2243   if (!cntrbox_pointin (c, p))
2244     return FALSE;
2245   info.f = 0;
2246   info.p[0] = ray.X1 = p[0];
2247   info.p[1] = ray.Y1 = p[1];
2248   ray.X2 = 0x7fffffff;
2249   ray.Y2 = p[1] + 1;
2250   if (setjmp (info.env) == 0)
2251     r_search (c->tree, &ray, NULL, crossing, &info);
2252   return info.f;
2253 }
2254
2255 BOOLp
2256 poly_CheckInside (POLYAREA * p, Vector v0)
2257 {
2258   PLINE *cur;
2259
2260   if ((p == NULL) || (v0 == NULL) || (p->contours == NULL))
2261     return FALSE;
2262   cur = p->contours;
2263   if (poly_InsideContour (cur, v0))
2264     {
2265       for (cur = cur->next; cur != NULL; cur = cur->next)
2266         if (poly_InsideContour (cur, v0))
2267           return FALSE;
2268       return TRUE;
2269     }
2270   return FALSE;
2271 }
2272
2273 BOOLp
2274 poly_M_CheckInside (POLYAREA * p, Vector v0)
2275 {
2276   POLYAREA *cur;
2277
2278   if ((p == NULL) || (v0 == NULL))
2279     return FALSE;
2280   cur = p;
2281   do
2282     {
2283       if (poly_CheckInside (cur, v0))
2284         return TRUE;
2285     }
2286   while ((cur = cur->f) != p);
2287   return FALSE;
2288 }
2289
2290 int
2291 poly_ContourInContour (PLINE * poly, PLINE * inner)
2292 {
2293   assert (poly != NULL);
2294   assert (inner != NULL);
2295   if (cntrbox_inside (inner, poly))
2296     return poly_InsideContour (poly, inner->head.point);
2297   return 0;
2298 }
2299
2300 void
2301 poly_Init (POLYAREA * p)
2302 {
2303   p->f = p->b = p;
2304   p->contours = NULL;
2305 }
2306
2307 POLYAREA *
2308 poly_Create (void)
2309 {
2310   POLYAREA *res;
2311
2312   if ((res = malloc (sizeof (POLYAREA))) != NULL)
2313     poly_Init (res);
2314   return res;
2315 }
2316
2317 void
2318 poly_FreeContours (PLINE **pline)
2319 {
2320   PLINE *pl;
2321
2322   while ((pl = *pline) != NULL)
2323     {
2324       *pline = pl->next;
2325       poly_DelContour (&pl);
2326     }
2327 }
2328
2329 void
2330 poly_Free (POLYAREA ** p)
2331 {
2332   POLYAREA *cur;
2333
2334   if (*p == NULL)
2335     return;
2336   for (cur = (*p)->f; cur != *p; cur = (*p)->f)
2337     {
2338       poly_FreeContours (&cur->contours);
2339       cur->f->b = cur->b;
2340       cur->b->f = cur->f;
2341       free (cur);
2342     }
2343   poly_FreeContours (&cur->contours);
2344   free (*p), *p = NULL;
2345 }
2346
2347 static BOOLp
2348 inside_sector (VNODE * pn, Vector p2)
2349 {
2350   Vector cdir, ndir, pdir;
2351   int p_c, n_c, p_n;
2352
2353   assert (pn != NULL);
2354   vect_sub (cdir, p2, pn->point);
2355   vect_sub (pdir, pn->point, pn->prev->point);
2356   vect_sub (ndir, pn->next->point, pn->point);
2357
2358   p_c = vect_det2 (pdir, cdir) >= 0;
2359   n_c = vect_det2 (ndir, cdir) >= 0;
2360   p_n = vect_det2 (pdir, ndir) >= 0;
2361
2362   if ((p_n && p_c && n_c) || ((!p_n) && (p_c || n_c)))
2363     return TRUE;
2364   else
2365     return FALSE;
2366 }                               /* inside_sector */
2367
2368 /* returns TRUE if bad contour */
2369 BOOLp
2370 poly_ChkContour (PLINE * a)
2371 {
2372   VNODE *a1, *a2, *hit1, *hit2;
2373   Vector i1, i2;
2374   int icnt;
2375
2376   assert (a != NULL);
2377   a1 = &a->head;
2378   do
2379     {
2380       a2 = a1;
2381       do
2382         {
2383           if (!node_neighbours (a1, a2) &&
2384               (icnt = vect_inters2 (a1->point, a1->next->point,
2385                                     a2->point, a2->next->point, i1, i2)) > 0)
2386             {
2387               if (icnt > 1)
2388                 return TRUE;
2389
2390               if (vect_dist2 (i1, a1->point) < EPSILON)
2391                 hit1 = a1;
2392               else if (vect_dist2 (i1, a1->next->point) < EPSILON)
2393                 hit1 = a1->next;
2394               else
2395                 return TRUE;
2396
2397               if (vect_dist2 (i1, a2->point) < EPSILON)
2398                 hit2 = a2;
2399               else if (vect_dist2 (i1, a2->next->point) < EPSILON)
2400                 hit2 = a2->next;
2401               else
2402                 return TRUE;
2403
2404 #if 1
2405               /* now check if they are inside each other */
2406               if (inside_sector (hit1, hit2->prev->point) ||
2407                   inside_sector (hit1, hit2->next->point) ||
2408                   inside_sector (hit2, hit1->prev->point) ||
2409                   inside_sector (hit2, hit1->next->point))
2410                 return TRUE;
2411 #endif
2412             }
2413         }
2414       while ((a2 = a2->next) != &a->head);
2415     }
2416   while ((a1 = a1->next) != &a->head);
2417   return FALSE;
2418 }
2419
2420
2421 BOOLp
2422 poly_Valid (POLYAREA * p)
2423 {
2424   PLINE *c;
2425
2426   if ((p == NULL) || (p->contours == NULL))
2427     return FALSE;
2428
2429   if (p->contours->Flags.orient == PLF_INV || poly_ChkContour (p->contours))
2430     {
2431 #ifndef NDEBUG
2432       VNODE *v;
2433       DEBUGP ("Invalid Outer PLINE\n");
2434       if (p->contours->Flags.orient == PLF_INV)
2435         DEBUGP ("failed orient\n");
2436       if (poly_ChkContour (p->contours))
2437         DEBUGP ("failed self-intersection\n");
2438       v = &p->contours->head;
2439       do
2440         {
2441           fprintf (stderr, "%d %d 100 100 \"\"]\n", v->point[0], v->point[1]);
2442           fprintf (stderr, "Line [%d %d ", v->point[0], v->point[1]);
2443         }
2444       while ((v = v->next) != &p->contours->head);
2445 #endif
2446       return FALSE;
2447     }
2448   for (c = p->contours->next; c != NULL; c = c->next)
2449     {
2450       if (c->Flags.orient == PLF_DIR ||
2451           poly_ChkContour (c) || !poly_ContourInContour (p->contours, c))
2452         {
2453 #ifndef NDEBUG
2454           VNODE *v;
2455           DEBUGP ("Invalid Inner PLINE orient = %d\n", c->Flags.orient);
2456           if (c->Flags.orient == PLF_DIR)
2457             DEBUGP ("failed orient\n");
2458           if (poly_ChkContour (c))
2459             DEBUGP ("failed self-intersection\n");
2460           if (!poly_ContourInContour (p->contours, c))
2461             DEBUGP ("failed containment\n");
2462           v = &c->head;
2463           do
2464             {
2465               fprintf (stderr, "%d %d 100 100 \"\"]\n", v->point[0],
2466                        v->point[1]);
2467               fprintf (stderr, "Line [%d %d ", v->point[0], v->point[1]);
2468             }
2469           while ((v = v->next) != &c->head);
2470 #endif
2471           return FALSE;
2472         }
2473     }
2474   return TRUE;
2475 }
2476
2477
2478 Vector vect_zero = { (long) 0, (long) 0 };
2479
2480 /*********************************************************************/
2481 /*             L o n g   V e c t o r   S t u f f                     */
2482 /*********************************************************************/
2483
2484 void
2485 vect_init (Vector v, double x, double y)
2486 {
2487   v[0] = (long) x;
2488   v[1] = (long) y;
2489 }                               /* vect_init */
2490
2491 #define Vzero(a)   ((a)[0] == 0. && (a)[1] == 0.)
2492
2493 #define Vsub(a,b,c) {(a)[0]=(b)[0]-(c)[0];(a)[1]=(b)[1]-(c)[1];}
2494
2495 int
2496 vect_equal (Vector v1, Vector v2)
2497 {
2498   return (v1[0] == v2[0] && v1[1] == v2[1]);
2499 }                               /* vect_equal */
2500
2501
2502 void
2503 vect_sub (Vector res, Vector v1, Vector v2)
2504 {
2505 Vsub (res, v1, v2)}             /* vect_sub */
2506
2507 void
2508 vect_min (Vector v1, Vector v2, Vector v3)
2509 {
2510   v1[0] = (v2[0] < v3[0]) ? v2[0] : v3[0];
2511   v1[1] = (v2[1] < v3[1]) ? v2[1] : v3[1];
2512 }                               /* vect_min */
2513
2514 void
2515 vect_max (Vector v1, Vector v2, Vector v3)
2516 {
2517   v1[0] = (v2[0] > v3[0]) ? v2[0] : v3[0];
2518   v1[1] = (v2[1] > v3[1]) ? v2[1] : v3[1];
2519 }                               /* vect_max */
2520
2521 double
2522 vect_len2 (Vector v)
2523 {
2524   return ((double) v[0] * v[0] + (double) v[1] * v[1]); /* why sqrt? only used for compares */
2525 }
2526
2527 double
2528 vect_dist2 (Vector v1, Vector v2)
2529 {
2530   double dx = v1[0] - v2[0];
2531   double dy = v1[1] - v2[1];
2532
2533   return (dx * dx + dy * dy);   /* why sqrt */
2534 }
2535
2536 /* value has sign of angle between vectors */
2537 double
2538 vect_det2 (Vector v1, Vector v2)
2539 {
2540   return (((double) v1[0] * v2[1]) - ((double) v2[0] * v1[1]));
2541 }
2542
2543 static double
2544 vect_m_dist (Vector v1, Vector v2)
2545 {
2546   double dx = v1[0] - v2[0];
2547   double dy = v1[1] - v2[1];
2548   double dd = (dx * dx + dy * dy);      /* sqrt */
2549
2550   if (dx > 0)
2551     return +dd;
2552   if (dx < 0)
2553     return -dd;
2554   if (dy > 0)
2555     return +dd;
2556   return -dd;
2557 }                               /* vect_m_dist */
2558
2559 /*
2560 vect_inters2
2561  (C) 1993 Klamer Schutte
2562  (C) 1997 Michael Leonov, Alexey Nikitin
2563 */
2564
2565 int
2566 vect_inters2 (Vector p1, Vector p2, Vector q1, Vector q2,
2567               Vector S1, Vector S2)
2568 {
2569   double s, t, deel;
2570   double rpx, rpy, rqx, rqy;
2571
2572   if (max (p1[0], p2[0]) < min (q1[0], q2[0]) ||
2573       max (q1[0], q2[0]) < min (p1[0], p2[0]) ||
2574       max (p1[1], p2[1]) < min (q1[1], q2[1]) ||
2575       max (q1[1], q2[1]) < min (p1[1], p2[1]))
2576     return 0;
2577
2578   rpx = p2[0] - p1[0];
2579   rpy = p2[1] - p1[1];
2580   rqx = q2[0] - q1[0];
2581   rqy = q2[1] - q1[1];
2582
2583   deel = rpy * rqx - rpx * rqy; /* -vect_det(rp,rq); */
2584
2585   /* coordinates are 30-bit integers so deel will be exactly zero
2586    * if the lines are parallel
2587    */
2588
2589   if (deel == 0)                /* parallel */
2590     {
2591       double dc1, dc2, d1, d2, h;       /* Check to see whether p1-p2 and q1-q2 are on the same line */
2592       Vector hp1, hq1, hp2, hq2, q1p1, q1q2;
2593
2594       Vsub2 (q1p1, q1, p1);
2595       Vsub2 (q1q2, q1, q2);
2596
2597
2598       /* If this product is not zero then p1-p2 and q1-q2 are not on same line! */
2599       if (vect_det2 (q1p1, q1q2) != 0)
2600         return 0;
2601       dc1 = 0;                  /* m_len(p1 - p1) */
2602
2603       dc2 = vect_m_dist (p1, p2);
2604       d1 = vect_m_dist (p1, q1);
2605       d2 = vect_m_dist (p1, q2);
2606
2607 /* Sorting the independent points from small to large */
2608       Vcpy2 (hp1, p1);
2609       Vcpy2 (hp2, p2);
2610       Vcpy2 (hq1, q1);
2611       Vcpy2 (hq2, q2);
2612       if (dc1 > dc2)
2613         {                       /* hv and h are used as help-variable. */
2614           Vswp2 (hp1, hp2);
2615           h = dc1, dc1 = dc2, dc2 = h;
2616         }
2617       if (d1 > d2)
2618         {
2619           Vswp2 (hq1, hq2);
2620           h = d1, d1 = d2, d2 = h;
2621         }
2622
2623 /* Now the line-pieces are compared */
2624
2625       if (dc1 < d1)
2626         {
2627           if (dc2 < d1)
2628             return 0;
2629           if (dc2 < d2)
2630             {
2631               Vcpy2 (S1, hp2);
2632               Vcpy2 (S2, hq1);
2633             }
2634           else
2635             {
2636               Vcpy2 (S1, hq1);
2637               Vcpy2 (S2, hq2);
2638             };
2639         }
2640       else
2641         {
2642           if (dc1 > d2)
2643             return 0;
2644           if (dc2 < d2)
2645             {
2646               Vcpy2 (S1, hp1);
2647               Vcpy2 (S2, hp2);
2648             }
2649           else
2650             {
2651               Vcpy2 (S1, hp1);
2652               Vcpy2 (S2, hq2);
2653             };
2654         }
2655       return (Vequ2 (S1, S2) ? 1 : 2);
2656     }
2657   else
2658     {                           /* not parallel */
2659       /*
2660        * We have the lines:
2661        * l1: p1 + s(p2 - p1)
2662        * l2: q1 + t(q2 - q1)
2663        * And we want to know the intersection point.
2664        * Calculate t:
2665        * p1 + s(p2-p1) = q1 + t(q2-q1)
2666        * which is similar to the two equations:
2667        * p1x + s * rpx = q1x + t * rqx
2668        * p1y + s * rpy = q1y + t * rqy
2669        * Multiplying these by rpy resp. rpx gives:
2670        * rpy * p1x + s * rpx * rpy = rpy * q1x + t * rpy * rqx
2671        * rpx * p1y + s * rpx * rpy = rpx * q1y + t * rpx * rqy
2672        * Subtracting these gives:
2673        * rpy * p1x - rpx * p1y = rpy * q1x - rpx * q1y + t * ( rpy * rqx - rpx * rqy )
2674        * So t can be isolated:
2675        * t = (rpy * ( p1x - q1x ) + rpx * ( - p1y + q1y )) / ( rpy * rqx - rpx * rqy )
2676        * and s can be found similarly
2677        * s = (rqy * (q1x - p1x) + rqx * (p1y - q1y))/( rqy * rpx - rqx * rpy)
2678        */
2679
2680       if (Vequ2 (q1, p1) || Vequ2 (q1, p2))
2681         {
2682           S1[0] = q1[0];
2683           S1[1] = q1[1];
2684         }
2685       else if (Vequ2 (q2, p1) || Vequ2 (q2, p2))
2686         {
2687           S1[0] = q2[0];
2688           S1[1] = q2[1];
2689         }
2690       else
2691         {
2692           s = (rqy * (p1[0] - q1[0]) + rqx * (q1[1] - p1[1])) / deel;
2693           if (s < 0 || s > 1.)
2694             return 0;
2695           t = (rpy * (p1[0] - q1[0]) + rpx * (q1[1] - p1[1])) / deel;
2696           if (t < 0 || t > 1.)
2697             return 0;
2698
2699           S1[0] = q1[0] + ROUND (t * rqx);
2700           S1[1] = q1[1] + ROUND (t * rqy);
2701         }
2702       return 1;
2703     }
2704 }                               /* vect_inters2 */
2705
2706 /* how about expanding polygons so that edges can be arcs rather than
2707  * lines. Consider using the third coordinate to store the radius of the
2708  * arc. The arc would pass through the vertex points. Positive radius
2709  * would indicate the arc bows left (center on right of P1-P2 path)
2710  * negative radius would put the center on the other side. 0 radius
2711  * would mean the edge is a line instead of arc
2712  * The intersections of the two circles centered at the vertex points
2713  * would determine the two possible arc centers. If P2.x > P1.x then
2714  * the center with smaller Y is selected for positive r. If P2.y > P1.y
2715  * then the center with greate X is selected for positive r.
2716  *
2717  * the vec_inters2() routine would then need to handle line-line
2718  * line-arc and arc-arc intersections.
2719  *
2720  * perhaps reverse tracing the arc would require look-ahead to check
2721  * for arcs
2722  */