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