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
  58 /*********************************************************************/
  59 /*              L o n g   V e c t o r   S t u f f                    */
  60 /*********************************************************************/
  61
  62 #define Vcopy(a,b) {(a)[0]=(b)[0];(a)[1]=(b)[1];}
  63 int vect_equal (Vector v1, Vector v2);
  64 void vect_init (Vector v, double x, double y);
  65 void vect_sub (Vector res, Vector v2, Vector v3);
  66
  67 void vect_min (Vector res, Vector v2, Vector v3);
  68 void vect_max (Vector res, Vector v2, Vector v3);
  69
  70 double vect_dist2 (Vector v1, Vector v2);
  71 double vect_det2 (Vector v1, Vector v2);
  72 double vect_len2 (Vector v1);
  73
  74 int vect_inters2 (Vector A, Vector B, Vector C, Vector D, Vector S1,
  75                   Vector S2);
  76
  77 /* note that a vertex v's Flags.status represents the edge defined by
  78  * v to v->next (i.e. the edge is forward of v)
  79  */
  80 #define ISECTED 3
  81 #define UNKNWN  0
  82 #define INSIDE  1
  83 #define OUTSIDE 2
  84 #define SHARED 3
  85 #define SHARED2 4
  86
  87 #define TOUCHES 99
  88
  89 #define NODE_LABEL(n)  ((n)->Flags.status)
  90 #define LABEL_NODE(n,l) ((n)->Flags.status = (l))
  91
  92 #define error(code)  longjmp(*(e), code)
  93
  94 #define MemGet(ptr, type) \
  95   if (UNLIKELY (((ptr) = (type *)malloc(sizeof(type))) == NULL))        \
  96     error(err_no_memory);
  97
  98 #undef DEBUG_LABEL
  99 #undef DEBUG_ALL_LABELS
 100 #undef DEBUG_JUMP
 101 #undef DEBUG_GATHER
 102 #undef DEBUG_ANGLE
 103 #undef DEBUG
 104 #ifdef DEBUG
 105 #define DEBUGP(...) pcb_fprintf(stderr, ## __VA_ARGS__)
 106 #else
 107 #define DEBUGP(...)
 108 #endif
 109
 110 /* ///////////////////////////////////////////////////////////////////////////// * /
 111 / *  2-Dimentional stuff
 112 / * ///////////////////////////////////////////////////////////////////////////// */
 113
 114 #define Vsub2(r,a,b)    {(r)[0] = (a)[0] - (b)[0]; (r)[1] = (a)[1] - (b)[1];}
 115 #define Vadd2(r,a,b)    {(r)[0] = (a)[0] + (b)[0]; (r)[1] = (a)[1] + (b)[1];}
 116 #define Vsca2(r,a,s)    {(r)[0] = (a)[0] * (s); (r)[1] = (a)[1] * (s);}
 117 #define Vcpy2(r,a)              {(r)[0] = (a)[0]; (r)[1] = (a)[1];}
 118 #define Vequ2(a,b)              ((a)[0] == (b)[0] && (a)[1] == (b)[1])
 119 #define Vadds(r,a,b,s)  {(r)[0] = ((a)[0] + (b)[0]) * (s); (r)[1] = ((a)[1] + (b)[1]) * (s);}
 120 #define Vswp2(a,b) { long t; \
 121         t = (a)[0], (a)[0] = (b)[0], (b)[0] = t; \
 122         t = (a)[1], (a)[1] = (b)[1], (b)[1] = t; \
 123 }
 124
 125 #ifdef DEBUG
 126 static char *theState (VNODE * v);
 127
 128 static void
 129 pline_dump (VNODE * v)
 130 {
 131   VNODE *s, *n;
 132
 133   s = v;
 134   do
 135     {
 136       n = v->next;
 137       pcb_fprintf (stderr, "Line [%#mS %#mS %#mS %#mS 10 10 \"%s\"]\n",
 138                v->point[0], v->point[1],
 139                n->point[0], n->point[1], theState (v));
 140     }
 141   while ((v = v->next) != s);
 142 }
 143
 144 static void
 145 poly_dump (POLYAREA * p)
 146 {
 147   POLYAREA *f = p;
 148   PLINE *pl;
 149
 150   do
 151     {
 152       pl = p->contours;
 153       do
 154         {
 155           pline_dump (&pl->head);
 156           fprintf (stderr, "NEXT PLINE\n");
 157         }
 158       while ((pl = pl->next) != NULL);
 159       fprintf (stderr, "NEXT POLY\n");
 160     }
 161   while ((p = p->f) != f);
 162 }
 163 #endif
 164
 165 /***************************************************************/
 166 /* routines for processing intersections */
 167
 168 /*
 169 node_add
 170  (C) 1993 Klamer Schutte
 171  (C) 1997 Alexey Nikitin, Michael Leonov
 172  (C) 2006 harry eaton
 173
 174  returns a bit field in new_point that indicates where the
 175  point was.
 176  1 means a new node was created and inserted
 177  4 means the intersection was not on the dest point
 178 */
 179 static VNODE *
 180 node_add_single (VNODE * dest, Vector po)
 181 {
 182   VNODE *p;
 183
 184   if (vect_equal (po, dest->point))
 185     return dest;
 186   if (vect_equal (po, dest->next->point))
 187     return dest->next;
 188   p = poly_CreateNode (po);
 189   if (p == NULL)
 190     return NULL;
 191   p->cvc_prev = p->cvc_next = NULL;
 192   p->Flags.status = UNKNWN;
 193   return p;
 194 }                               /* node_add */
 195
 196 #define ISECT_BAD_PARAM (-1)
 197 #define ISECT_NO_MEMORY (-2)
 198
 199
 200 /*
 201 new_descriptor
 202   (C) 2006 harry eaton
 203 */
 204 static CVCList *
 205 new_descriptor (VNODE * a, char poly, char side)
 206 {
 207   CVCList *l = (CVCList *) malloc (sizeof (CVCList));
 208   Vector v;
 209   register double ang, dx, dy;
 210
 211   if (!l)
 212     return NULL;
 213   l->head = NULL;
 214   l->parent = a;
 215   l->poly = poly;
 216   l->side = side;
 217   l->next = l->prev = l;
 218   if (side == 'P')              /* previous */
 219     vect_sub (v, a->prev->point, a->point);
 220   else                          /* next */
 221     vect_sub (v, a->next->point, a->point);
 222   /* Uses slope/(slope+1) in quadrant 1 as a proxy for the angle.
 223    * It still has the same monotonic sort result
 224    * and is far less expensive to compute than the real angle.
 225    */
 226   if (vect_equal (v, vect_zero))
 227     {
 228       if (side == 'P')
 229         {
 230           if (a->prev->cvc_prev == (CVCList *) - 1)
 231             a->prev->cvc_prev = a->prev->cvc_next = NULL;
 232           poly_ExclVertex (a->prev);
 233           vect_sub (v, a->prev->point, a->point);
 234         }
 235       else
 236         {
 237           if (a->next->cvc_prev == (CVCList *) - 1)
 238             a->next->cvc_prev = a->next->cvc_next = NULL;
 239           poly_ExclVertex (a->next);
 240           vect_sub (v, a->next->point, a->point);
 241         }
 242     }
 243   assert (!vect_equal (v, vect_zero));
 244   dx = fabs ((double) v[0]);
 245   dy = fabs ((double) v[1]);
 246   ang = dy / (dy + dx);
 247   /* now move to the actual quadrant */
 248   if (v[0] < 0 && v[1] >= 0)
 249     ang = 2.0 - ang;            /* 2nd quadrant */
 250   else if (v[0] < 0 && v[1] < 0)
 251     ang += 2.0;                 /* 3rd quadrant */
 252   else if (v[0] >= 0 && v[1] < 0)
 253     ang = 4.0 - ang;            /* 4th quadrant */
 254   l->angle = ang;
 255   assert (ang >= 0.0 && ang <= 4.0);
 256 #ifdef DEBUG_ANGLE
 257   DEBUGP ("node on %c at %#mD assigned angle %g on side %c\n", poly,
 258           a->point[0], a->point[1], ang, side);
 259 #endif
 260   return l;
 261 }
 262
 263 /*
 264 insert_descriptor
 265   (C) 2006 harry eaton
 266
 267    argument a is a cross-vertex node.
 268    argument poly is the polygon it comes from ('A' or 'B')
 269    argument side is the side this descriptor goes on ('P' for previous
 270    'N' for next.
 271    argument start is the head of the list of cvclists
 272 */
 273 static CVCList *
 274 insert_descriptor (VNODE * a, char poly, char side, CVCList * start)
 275 {
 276   CVCList *l, *newone, *big, *small;
 277
 278   if (!(newone = new_descriptor (a, poly, side)))
 279     return NULL;
 280   /* search for the CVCList for this point */
 281   if (!start)
 282     {
 283       start = newone;           /* return is also new, so we know where start is */
 284       start->head = newone;     /* circular list */
 285       return newone;
 286     }
 287   else
 288     {
 289       l = start;
 290       do
 291         {
 292           assert (l->head);
 293           if (l->parent->point[0] == a->point[0]
 294               && l->parent->point[1] == a->point[1])
 295             {                   /* this CVCList is at our point */
 296               start = l;
 297               newone->head = l->head;
 298               break;
 299             }
 300           if (l->head->parent->point[0] == start->parent->point[0]
 301               && l->head->parent->point[1] == start->parent->point[1])
 302             {
 303               /* this seems to be a new point */
 304               /* link this cvclist to the list of all cvclists */
 305               for (; l->head != newone; l = l->next)
 306                 l->head = newone;
 307               newone->head = start;
 308               return newone;
 309             }
 310           l = l->head;
 311         }
 312       while (1);
 313     }
 314   assert (start);
 315   l = big = small = start;
 316   do
 317     {
 318       if (l->next->angle < l->angle)    /* find start/end of list */
 319         {
 320           small = l->next;
 321           big = l;
 322         }
 323       else if (newone->angle >= l->angle && newone->angle <= l->next->angle)
 324         {
 325           /* insert new cvc if it lies between existing points */
 326           newone->prev = l;
 327           newone->next = l->next;
 328           l->next = l->next->prev = newone;
 329           return newone;
 330         }
 331     }
 332   while ((l = l->next) != start);
 333   /* didn't find it between points, it must go on an end */
 334   if (big->angle <= newone->angle)
 335     {
 336       newone->prev = big;
 337       newone->next = big->next;
 338       big->next = big->next->prev = newone;
 339       return newone;
 340     }
 341   assert (small->angle >= newone->angle);
 342   newone->next = small;
 343   newone->prev = small->prev;
 344   small->prev = small->prev->next = newone;
 345   return newone;
 346 }
 347
 348 /*
 349 node_add_point
 350  (C) 1993 Klamer Schutte
 351  (C) 1997 Alexey Nikitin, Michael Leonov
 352
 353  return 1 if new node in b, 2 if new node in a and 3 if new node in both
 354 */
 355
 356 static VNODE *
 357 node_add_single_point (VNODE * a, Vector p)
 358 {
 359   VNODE *next_a, *new_node;
 360
 361   next_a = a->next;
 362
 363   new_node = node_add_single (a, p);
 364   assert (new_node != NULL);
 365
 366   new_node->cvc_prev = new_node->cvc_next = (CVCList *) - 1;
 367
 368   if (new_node == a || new_node == next_a)
 369     return NULL;
 370
 371   return new_node;
 372 }                               /* node_add_point */
 373
 374 /*
 375 node_label
 376  (C) 2006 harry eaton
 377 */
 378 static unsigned int
 379 node_label (VNODE * pn)
 380 {
 381   CVCList *first_l, *l;
 382   char this_poly;
 383   int region = UNKNWN;
 384
 385   assert (pn);
 386   assert (pn->cvc_next);
 387   this_poly = pn->cvc_next->poly;
 388   /* search counter-clockwise in the cross vertex connectivity (CVC) list
 389    *
 390    * check for shared edges (that could be prev or next in the list since the angles are equal)
 391    * and check if this edge (pn -> pn->next) is found between the other poly's entry and exit
 392    */
 393   if (pn->cvc_next->angle == pn->cvc_next->prev->angle)
 394     l = pn->cvc_next->prev;
 395   else
 396     l = pn->cvc_next->next;
 397
 398   first_l = l;
 399   while ((l->poly == this_poly) && (l != first_l->prev))
 400     l = l->next;
 401   assert (l->poly != this_poly);
 402
 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   int intersected;
 496 } seg;
 497
 498 typedef struct _insert_node_task insert_node_task;
 499
 500 struct _insert_node_task
 501 {
 502   insert_node_task *next;
 503   seg * node_seg;
 504   VNODE *new_node;
 505 };
 506
 507 typedef struct info
 508 {
 509   double m, b;
 510   rtree_t *tree;
 511   VNODE *v;
 512   struct seg *s;
 513   jmp_buf *env, sego, *touch;
 514   int need_restart;
 515   insert_node_task *node_insert_list;
 516 } info;
 517
 518 typedef struct contour_info
 519 {
 520   PLINE *pa;
 521   jmp_buf restart;
 522   jmp_buf *getout;
 523   int need_restart;
 524   insert_node_task *node_insert_list;
 525 } contour_info;
 526
 527
 528 /*
 529  * adjust_tree()
 530  * (C) 2006 harry eaton
 531  * This replaces the segment in the tree with the two new segments after
 532  * a vertex has been added
 533  */
 534 static int
 535 adjust_tree (rtree_t * tree, struct seg *s)
 536 {
 537   struct seg *q;
 538
 539   q = (seg *)malloc (sizeof (struct seg));
 540   if (!q)
 541     return 1;
 542   q->intersected = 0;
 543   q->v = s->v;
 544   q->p = s->p;
 545   q->box.X1 = min (q->v->point[0], q->v->next->point[0]);
 546   q->box.X2 = max (q->v->point[0], q->v->next->point[0]) + 1;
 547   q->box.Y1 = min (q->v->point[1], q->v->next->point[1]);
 548   q->box.Y2 = max (q->v->point[1], q->v->next->point[1]) + 1;
 549   r_insert_entry (tree, (const BoxType *) q, 1);
 550   q = (seg *)malloc (sizeof (struct seg));
 551   if (!q)
 552     return 1;
 553   q->intersected = 0;
 554   q->v = s->v->next;
 555   q->p = s->p;
 556   q->box.X1 = min (q->v->point[0], q->v->next->point[0]);
 557   q->box.X2 = max (q->v->point[0], q->v->next->point[0]) + 1;
 558   q->box.Y1 = min (q->v->point[1], q->v->next->point[1]);
 559   q->box.Y2 = max (q->v->point[1], q->v->next->point[1]) + 1;
 560   r_insert_entry (tree, (const BoxType *) q, 1);
 561   r_delete_entry (tree, (const BoxType *) s);
 562   return 0;
 563 }
 564
 565 /*
 566  * seg_in_region()
 567  * (C) 2006, harry eaton
 568  * This prunes the search for boxes that don't intersect the segment.
 569  */
 570 static int
 571 seg_in_region (const BoxType * b, void *cl)
 572 {
 573   struct info *i = (struct info *) cl;
 574   double y1, y2;
 575   /* for zero slope the search is aligned on the axis so it is already pruned */
 576   if (i->m == 0.)
 577     return 1;
 578   y1 = i->m * b->X1 + i->b;
 579   y2 = i->m * b->X2 + i->b;
 580   if (min (y1, y2) >= b->Y2)
 581     return 0;
 582   if (max (y1, y2) < b->Y1)
 583     return 0;
 584   return 1;                     /* might intersect */
 585 }
 586
 587 /* Prepend a deferred node-insersion task to a list */
 588 static insert_node_task *
 589 prepend_insert_node_task (insert_node_task *list, seg *seg, VNODE *new_node)
 590 {
 591   insert_node_task *task = (insert_node_task *)malloc (sizeof (*task));
 592   task->node_seg = seg;
 593   task->new_node = new_node;
 594   task->next = list;
 595   return task;
 596 }
 597
 598 /*
 599  * seg_in_seg()
 600  * (C) 2006 harry eaton
 601  * This routine checks if the segment in the tree intersect the search segment.
 602  * If it does, the plines are marked as intersected and the point is marked for
 603  * the cvclist. If the point is not already a vertex, a new vertex is inserted
 604  * and the search for intersections starts over at the beginning.
 605  * That is potentially a significant time penalty, but it does solve the snap rounding
 606  * problem. There are efficient algorithms for finding intersections with snap
 607  * rounding, but I don't have time to implement them right now.
 608  */
 609 static int
 610 seg_in_seg (const BoxType * b, void *cl)
 611 {
 612   struct info *i = (struct info *) cl;
 613   struct seg *s = (struct seg *) b;
 614   Vector s1, s2;
 615   int cnt;
 616   VNODE *new_node;
 617
 618   /* When new nodes are added at the end of a pass due to an intersection
 619    * the segments may be altered. If either segment we're looking at has
 620    * already been intersected this pass, skip it until the next pass.
 621    */
 622   if (s->intersected || i->s->intersected)
 623     return 0;
 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       bool done_insert_on_i = false;
 636       new_node = node_add_single_point (i->v, cnt > 1 ? s2 : s1);
 637       if (new_node != NULL)
 638         {
 639 #ifdef DEBUG_INTERSECT
 640           DEBUGP ("new intersection on segment \"i\" at %#mD\n",
 641                   cnt > 1 ? s2[0] : s1[0], cnt > 1 ? s2[1] : s1[1]);
 642 #endif
 643           i->node_insert_list =
 644             prepend_insert_node_task (i->node_insert_list, i->s, new_node);
 645           i->s->intersected = 1;
 646           done_insert_on_i = true;
 647         }
 648       new_node = node_add_single_point (s->v, cnt > 1 ? s2 : s1);
 649       if (new_node != NULL)
 650         {
 651 #ifdef DEBUG_INTERSECT
 652           DEBUGP ("new intersection on segment \"s\" at %#mD\n",
 653                   cnt > 1 ? s2[0] : s1[0], cnt > 1 ? s2[1] : s1[1]);
 654 #endif
 655           i->node_insert_list =
 656             prepend_insert_node_task (i->node_insert_list, s, new_node);
 657           s->intersected = 1;
 658           return 0; /* Keep looking for intersections with segment "i" */
 659         }
 660       /* Skip any remaining r_search hits against segment i, as any futher
 661        * intersections will be rejected until the next pass anyway.
 662        */
 663       if (done_insert_on_i)
 664         longjmp (*i->env, 1);
 665     }
 666   return 0;
 667 }
 668
 669 static void *
 670 make_edge_tree (PLINE * pb)
 671 {
 672   struct seg *s;
 673   VNODE *bv;
 674   rtree_t *ans = r_create_tree (NULL, 0, 0);
 675   bv = &pb->head;
 676   do
 677     {
 678       s = (seg *)malloc (sizeof (struct seg));
 679       s->intersected = 0;
 680       if (bv->point[0] < bv->next->point[0])
 681         {
 682           s->box.X1 = bv->point[0];
 683           s->box.X2 = bv->next->point[0] + 1;
 684         }
 685       else
 686         {
 687           s->box.X2 = bv->point[0] + 1;
 688           s->box.X1 = bv->next->point[0];
 689         }
 690       if (bv->point[1] < bv->next->point[1])
 691         {
 692           s->box.Y1 = bv->point[1];
 693           s->box.Y2 = bv->next->point[1] + 1;
 694         }
 695       else
 696         {
 697           s->box.Y2 = bv->point[1] + 1;
 698           s->box.Y1 = bv->next->point[1];
 699         }
 700       s->v = bv;
 701       s->p = pb;
 702       r_insert_entry (ans, (const BoxType *) s, 1);
 703     }
 704   while ((bv = bv->next) != &pb->head);
 705   return (void *) ans;
 706 }
 707
 708 static int
 709 get_seg (const BoxType * b, void *cl)
 710 {
 711   struct info *i = (struct info *) cl;
 712   struct seg *s = (struct seg *) b;
 713   if (i->v == s->v)
 714     {
 715       i->s = s;
 716       longjmp (i->sego, 1);
 717     }
 718   return 0;
 719 }
 720
 721 /*
 722  * intersect() (and helpers)
 723  * (C) 2006, harry eaton
 724  * This uses an rtree to find A-B intersections. Whenever a new vertex is
 725  * added, the search for intersections is re-started because the rounding
 726  * could alter the topology otherwise.
 727  * This should use a faster algorithm for snap rounding intersection finding.
 728  * The best algorthim is probably found in:
 729  *
 730  * "Improved output-sensitive snap rounding," John Hershberger, Proceedings
 731  * of the 22nd annual symposium on Computational geomerty, 2006, pp 357-366.
 732  * http://doi.acm.org/10.1145/1137856.1137909
 733  *
 734  * Algorithms described by de Berg, or Goodrich or Halperin, or Hobby would
 735  * probably work as well.
 736  *
 737  */
 738
 739 static int
 740 contour_bounds_touch (const BoxType * b, void *cl)
 741 {
 742   contour_info *c_info = (contour_info *) cl;
 743   PLINE *pa = c_info->pa;
 744   PLINE *pb = (PLINE *) b;
 745   PLINE *rtree_over;
 746   PLINE *looping_over;
 747   VNODE *av;                    /* node iterators */
 748   struct info info;
 749   BoxType box;
 750   jmp_buf restart;
 751
 752   /* Have seg_in_seg return to our desired location if it touches */
 753   info.env = &restart;
 754   info.touch = c_info->getout;
 755   info.need_restart = 0;
 756   info.node_insert_list = c_info->node_insert_list;
 757
 758   /* Pick which contour has the fewer points, and do the loop
 759    * over that. The r_tree makes hit-testing against a contour
 760    * faster, so we want to do that on the bigger contour.
 761    */
 762   if (pa->Count < pb->Count)
 763     {
 764       rtree_over = pb;
 765       looping_over = pa;
 766     }
 767   else
 768     {
 769       rtree_over = pa;
 770       looping_over = pb;
 771     }
 772
 773   av = &looping_over->head;
 774   do                            /* Loop over the nodes in the smaller contour */
 775     {
 776       /* check this edge for any insertions */
 777       double dx;
 778       info.v = av;
 779       /* compute the slant for region trimming */
 780       dx = av->next->point[0] - av->point[0];
 781       if (dx == 0)
 782         info.m = 0;
 783       else
 784         {
 785           info.m = (av->next->point[1] - av->point[1]) / dx;
 786           info.b = av->point[1] - info.m * av->point[0];
 787         }
 788       box.X2 = (box.X1 = av->point[0]) + 1;
 789       box.Y2 = (box.Y1 = av->point[1]) + 1;
 790
 791       /* fill in the segment in info corresponding to this node */
 792       if (setjmp (info.sego) == 0)
 793         {
 794           r_search (looping_over->tree, &box, NULL, get_seg, &info);
 795           assert (0);
 796         }
 797
 798       /* If we're going to have another pass anyway, skip this */
 799       if (info.s->intersected && info.node_insert_list != NULL)
 800         continue;
 801
 802       if (setjmp (restart))
 803         continue;
 804
 805       /* NB: If this actually hits anything, we are teleported back to the beginning */
 806       info.tree = rtree_over->tree;
 807       if (info.tree)
 808         if (UNLIKELY (r_search (info.tree, &info.s->box,
 809                                 seg_in_region, seg_in_seg, &info)))
 810           assert (0); /* XXX: Memory allocation failure */
 811     }
 812   while ((av = av->next) != &looping_over->head);
 813
 814   c_info->node_insert_list = info.node_insert_list;
 815   if (info.need_restart)
 816     c_info->need_restart = 1;
 817   return 0;
 818 }
 819
 820 static int
 821 intersect_impl (jmp_buf * jb, POLYAREA * b, POLYAREA * a, int add)
 822 {
 823   POLYAREA *t;
 824   PLINE *pa;
 825   contour_info c_info;
 826   int need_restart = 0;
 827   insert_node_task *task;
 828   c_info.need_restart = 0;
 829   c_info.node_insert_list = NULL;
 830
 831   /* Search the r-tree of the object with most contours
 832    * We loop over the contours of "a". Swap if necessary.
 833    */
 834   if (a->contour_tree->size > b->contour_tree->size)
 835     {
 836       t = b;
 837       b = a;
 838       a = t;
 839     }
 840
 841   for (pa = a->contours; pa; pa = pa->next)     /* Loop over the contours of POLYAREA "a" */
 842     {
 843       BoxType sb;
 844       jmp_buf out;
 845       int retval;
 846
 847       c_info.getout = NULL;
 848       c_info.pa = pa;
 849
 850       if (!add)
 851         {
 852           retval = setjmp (out);
 853           if (retval)
 854             {
 855               /* The intersection test short-circuited back here,
 856                * we need to clean up, then longjmp to jb */
 857               longjmp (*jb, retval);
 858             }
 859           c_info.getout = &out;
 860         }
 861
 862       sb.X1 = pa->xmin;
 863       sb.Y1 = pa->ymin;
 864       sb.X2 = pa->xmax + 1;
 865       sb.Y2 = pa->ymax + 1;
 866
 867       r_search (b->contour_tree, &sb, NULL, contour_bounds_touch, &c_info);
 868       if (c_info.need_restart)
 869         need_restart = 1;
 870     }
 871
 872   /* Process any deferred node insersions */
 873   task = c_info.node_insert_list;
 874   while (task != NULL)
 875     {
 876       insert_node_task *next = task->next;
 877
 878       /* Do insersion */
 879       task->new_node->prev = task->node_seg->v;
 880       task->new_node->next = task->node_seg->v->next;
 881       task->node_seg->v->next->prev = task->new_node;
 882       task->node_seg->v->next = task->new_node;
 883       task->node_seg->p->Count++;
 884
 885       cntrbox_adjust (task->node_seg->p, task->new_node->point);
 886       if (adjust_tree (task->node_seg->p->tree, task->node_seg))
 887         assert (0); /* XXX: Memory allocation failure */
 888
 889       need_restart = 1; /* Any new nodes could intersect */
 890
 891       free (task);
 892       task = next;
 893     }
 894
 895   return need_restart;
 896 }
 897
 898 static int
 899 intersect (jmp_buf * jb, POLYAREA * b, POLYAREA * a, int add)
 900 {
 901   int call_count = 1;
 902   while (intersect_impl (jb, b, a, add))
 903     call_count++;
 904   return 0;
 905 }
 906
 907 static void
 908 M_POLYAREA_intersect (jmp_buf * e, POLYAREA * afst, POLYAREA * bfst, int add)
 909 {
 910   POLYAREA *a = afst, *b = bfst;
 911   PLINE *curcA, *curcB;
 912   CVCList *the_list = NULL;
 913
 914   if (a == NULL || b == NULL)
 915     error (err_bad_parm);
 916   do
 917     {
 918       do
 919         {
 920           if (a->contours->xmax >= b->contours->xmin &&
 921               a->contours->ymax >= b->contours->ymin &&
 922               a->contours->xmin <= b->contours->xmax &&
 923               a->contours->ymin <= b->contours->ymax)
 924             {
 925               if (UNLIKELY (intersect (e, a, b, add)))
 926                 error (err_no_memory);
 927             }
 928         }
 929       while (add && (a = a->f) != afst);
 930       for (curcB = b->contours; curcB != NULL; curcB = curcB->next)
 931         if (curcB->Flags.status == ISECTED)
 932           {
 933             the_list = add_descriptors (curcB, 'B', the_list);
 934             if (UNLIKELY (the_list == NULL))
 935               error (err_no_memory);
 936           }
 937     }
 938   while (add && (b = b->f) != bfst);
 939   do
 940     {
 941       for (curcA = a->contours; curcA != NULL; curcA = curcA->next)
 942         if (curcA->Flags.status == ISECTED)
 943           {
 944             the_list = add_descriptors (curcA, 'A', the_list);
 945             if (UNLIKELY (the_list == NULL))
 946               error (err_no_memory);
 947           }
 948     }
 949   while (add && (a = a->f) != afst);
 950 }                               /* M_POLYAREA_intersect */
 951
 952 static inline int
 953 cntrbox_inside (PLINE * c1, PLINE * c2)
 954 {
 955   assert (c1 != NULL && c2 != NULL);
 956   return ((c1->xmin >= c2->xmin) &&
 957           (c1->ymin >= c2->ymin) &&
 958           (c1->xmax <= c2->xmax) && (c1->ymax <= c2->ymax));
 959 }
 960
 961 /*****************************************************************/
 962 /* Routines for making labels */
 963
 964 static int
 965 count_contours_i_am_inside (const BoxType * b, void *cl)
 966 {
 967   PLINE *me = (PLINE *) cl;
 968   PLINE *check = (PLINE *) b;
 969
 970   if (poly_ContourInContour (check, me))
 971     return 1;
 972   return 0;
 973 }
 974
 975 /* cntr_in_M_POLYAREA
 976 returns poly is inside outfst ? TRUE : FALSE */
 977 static int
 978 cntr_in_M_POLYAREA (PLINE * poly, POLYAREA * outfst, BOOLp test)
 979 {
 980   POLYAREA *outer = outfst;
 981   heap_t *heap;
 982
 983   assert (poly != NULL);
 984   assert (outer != NULL);
 985
 986   heap = heap_create ();
 987   do
 988     {
 989       if (cntrbox_inside (poly, outer->contours))
 990         heap_insert (heap, outer->contours->area, (void *) outer);
 991     }
 992   /* if checking touching, use only the first polygon */
 993   while (!test && (outer = outer->f) != outfst);
 994   /* we need only check the smallest poly container
 995    * but we must loop in case the box containter is not
 996    * the poly container */
 997   do
 998     {
 999       if (heap_is_empty (heap))
1000         break;
1001       outer = (POLYAREA *) heap_remove_smallest (heap);
1002
1003       switch (r_search
1004               (outer->contour_tree, (BoxType *) poly, NULL,
1005                count_contours_i_am_inside, poly))
1006         {
1007         case 0:         /* Didn't find anything in this piece, Keep looking */
1008           break;
1009         case 1:         /* Found we are inside this piece, and not any of its holes */
1010           heap_destroy (&heap);
1011           return TRUE;
1012         case 2:         /* Found inside a hole in the smallest polygon so far. No need to check the other polygons */
1013           heap_destroy (&heap);
1014           return FALSE;
1015         default:
1016           printf ("Something strange here\n");
1017           break;
1018         }
1019     }
1020   while (1);
1021   heap_destroy (&heap);
1022   return FALSE;
1023 }                               /* cntr_in_M_POLYAREA */
1024
1025 #ifdef DEBUG
1026
1027 static char *
1028 theState (VNODE * v)
1029 {
1030   static char u[] = "UNKNOWN";
1031   static char i[] = "INSIDE";
1032   static char o[] = "OUTSIDE";
1033   static char s[] = "SHARED";
1034   static char s2[] = "SHARED2";
1035
1036   switch (NODE_LABEL (v))
1037     {
1038     case INSIDE:
1039       return i;
1040     case OUTSIDE:
1041       return o;
1042     case SHARED:
1043       return s;
1044     case SHARED2:
1045       return s2;
1046     default:
1047       return u;
1048     }
1049 }
1050
1051 #ifdef DEBUG_ALL_LABELS
1052 static void
1053 print_labels (PLINE * a)
1054 {
1055   VNODE *c = &a->head;
1056
1057   do
1058     {
1059       DEBUGP ("%#mD->%#mD labeled %s\n", c->point[0], c->point[1],
1060               c->next->point[0], c->next->point[1], theState (c));
1061     }
1062   while ((c = c->next) != &a->head);
1063 }
1064 #endif
1065 #endif
1066
1067 /*
1068 label_contour
1069  (C) 2006 harry eaton
1070  (C) 1993 Klamer Schutte
1071  (C) 1997 Alexey Nikitin, Michael Leonov
1072 */
1073
1074 static BOOLp
1075 label_contour (PLINE * a)
1076 {
1077   VNODE *cur = &a->head;
1078   VNODE *first_labelled = NULL;
1079   int label = UNKNWN;
1080
1081   do
1082     {
1083       if (cur->cvc_next)        /* examine cross vertex */
1084         {
1085           label = node_label (cur);
1086           if (first_labelled == NULL)
1087             first_labelled = cur;
1088           continue;
1089         }
1090
1091       if (first_labelled == NULL)
1092         continue;
1093
1094       /* This labels nodes which aren't cross-connected */
1095       assert (label == INSIDE || label == OUTSIDE);
1096       LABEL_NODE (cur, label);
1097     }
1098   while ((cur = cur->next) != first_labelled);
1099 #ifdef DEBUG_ALL_LABELS
1100   print_labels (a);
1101   DEBUGP ("\n\n");
1102 #endif
1103   return FALSE;
1104 }                               /* label_contour */
1105
1106 static BOOLp
1107 cntr_label_POLYAREA (PLINE * poly, POLYAREA * ppl, BOOLp test)
1108 {
1109   assert (ppl != NULL && ppl->contours != NULL);
1110   if (poly->Flags.status == ISECTED)
1111     {
1112       label_contour (poly);     /* should never get here when BOOLp is true */
1113     }
1114   else if (cntr_in_M_POLYAREA (poly, ppl, test))
1115     {
1116       if (test)
1117         return TRUE;
1118       poly->Flags.status = INSIDE;
1119     }
1120   else
1121     {
1122       if (test)
1123         return false;
1124       poly->Flags.status = OUTSIDE;
1125     }
1126   return FALSE;
1127 }                               /* cntr_label_POLYAREA */
1128
1129 static BOOLp
1130 M_POLYAREA_label_separated (PLINE * afst, POLYAREA * b, BOOLp touch)
1131 {
1132   PLINE *curc = afst;
1133
1134   for (curc = afst; curc != NULL; curc = curc->next)
1135     {
1136       if (cntr_label_POLYAREA (curc, b, touch) && touch)
1137         return TRUE;
1138     }
1139   return FALSE;
1140 }
1141
1142 static BOOLp
1143 M_POLYAREA_label (POLYAREA * afst, POLYAREA * b, BOOLp touch)
1144 {
1145   POLYAREA *a = afst;
1146   PLINE *curc;
1147
1148   assert (a != NULL);
1149   do
1150     {
1151       for (curc = a->contours; curc != NULL; curc = curc->next)
1152         if (cntr_label_POLYAREA (curc, b, touch))
1153           {
1154             if (touch)
1155               return TRUE;
1156           }
1157     }
1158   while (!touch && (a = a->f) != afst);
1159   return FALSE;
1160 }
1161
1162 /****************************************************************/
1163
1164 /* routines for temporary storing resulting contours */
1165 static void
1166 InsCntr (jmp_buf * e, PLINE * c, POLYAREA ** dst)
1167 {
1168   POLYAREA *newp;
1169
1170   if (*dst == NULL)
1171     {
1172       MemGet (*dst, POLYAREA);
1173       (*dst)->f = (*dst)->b = *dst;
1174       newp = *dst;
1175     }
1176   else
1177     {
1178       MemGet (newp, POLYAREA);
1179       newp->f = *dst;
1180       newp->b = (*dst)->b;
1181       newp->f->b = newp->b->f = newp;
1182     }
1183   newp->contours = c;
1184   newp->contour_tree = r_create_tree (NULL, 0, 0);
1185   r_insert_entry (newp->contour_tree, (BoxType *) c, 0);
1186   c->next = NULL;
1187 }                               /* InsCntr */
1188
1189 static void
1190 PutContour (jmp_buf * e, PLINE * cntr, POLYAREA ** contours, PLINE ** holes,
1191             POLYAREA * owner, POLYAREA * parent, PLINE * parent_contour)
1192 {
1193   assert (cntr != NULL);
1194   assert (cntr->Count > 2);
1195   cntr->next = NULL;
1196
1197   if (cntr->Flags.orient == PLF_DIR)
1198     {
1199       if (owner != NULL)
1200         r_delete_entry (owner->contour_tree, (BoxType *) cntr);
1201       InsCntr (e, cntr, contours);
1202     }
1203   /* put hole into temporary list */
1204   else
1205     {
1206       /* if we know this belongs inside the parent, put it there now */
1207       if (parent_contour)
1208         {
1209           cntr->next = parent_contour->next;
1210           parent_contour->next = cntr;
1211           if (owner != parent)
1212             {
1213               if (owner != NULL)
1214                 r_delete_entry (owner->contour_tree, (BoxType *) cntr);
1215               r_insert_entry (parent->contour_tree, (BoxType *) cntr, 0);
1216             }
1217         }
1218       else
1219         {
1220           cntr->next = *holes;
1221           *holes = cntr;        /* let cntr be 1st hole in list */
1222           /* We don't insert the holes into an r-tree,
1223            * they just form a linked list */
1224           if (owner != NULL)
1225             r_delete_entry (owner->contour_tree, (BoxType *) cntr);
1226         }
1227     }
1228 }                               /* PutContour */
1229
1230 static inline void
1231 remove_contour (POLYAREA * piece, PLINE * prev_contour, PLINE * contour,
1232                 int remove_rtree_entry)
1233 {
1234   if (piece->contours == contour)
1235     piece->contours = contour->next;
1236   else if (prev_contour != NULL)
1237     {
1238       assert (prev_contour->next == contour);
1239       prev_contour->next = contour->next;
1240     }
1241
1242   contour->next = NULL;
1243
1244   if (remove_rtree_entry)
1245     r_delete_entry (piece->contour_tree, (BoxType *) contour);
1246 }
1247
1248 struct polyarea_info
1249 {
1250   BoxType BoundingBox;
1251   POLYAREA *pa;
1252 };
1253
1254 static int
1255 heap_it (const BoxType * b, void *cl)
1256 {
1257   heap_t *heap = (heap_t *) cl;
1258   struct polyarea_info *pa_info = (struct polyarea_info *) b;
1259   PLINE *p = pa_info->pa->contours;
1260   if (p->Count == 0)
1261     return 0;                   /* how did this happen? */
1262   heap_insert (heap, p->area, pa_info);
1263   return 1;
1264 }
1265
1266 struct find_inside_info
1267 {
1268   jmp_buf jb;
1269   PLINE *want_inside;
1270   PLINE *result;
1271 };
1272
1273 static int
1274 find_inside (const BoxType * b, void *cl)
1275 {
1276   struct find_inside_info *info = (struct find_inside_info *) cl;
1277   PLINE *check = (PLINE *) b;
1278   /* Do test on check to see if it inside info->want_inside */
1279   /* If it is: */
1280   if (check->Flags.orient == PLF_DIR)
1281     {
1282       return 0;
1283     }
1284   if (poly_ContourInContour (info->want_inside, check))
1285     {
1286       info->result = check;
1287       longjmp (info->jb, 1);
1288     }
1289   return 0;
1290 }
1291
1292 static void
1293 InsertHoles (jmp_buf * e, POLYAREA * dest, PLINE ** src)
1294 {
1295   POLYAREA *curc;
1296   PLINE *curh, *container;
1297   heap_t *heap;
1298   rtree_t *tree;
1299   int i;
1300   int num_polyareas = 0;
1301   struct polyarea_info *all_pa_info, *pa_info;
1302
1303   if (*src == NULL)
1304     return;                     /* empty hole list */
1305   if (dest == NULL)
1306     error (err_bad_parm);       /* empty contour list */
1307
1308   /* Count dest polyareas */
1309   curc = dest;
1310   do
1311     {
1312       num_polyareas++;
1313     }
1314   while ((curc = curc->f) != dest);
1315
1316   /* make a polyarea info table */
1317   /* make an rtree of polyarea info table */
1318   all_pa_info = (struct polyarea_info *) malloc (sizeof (struct polyarea_info) * num_polyareas);
1319   tree = r_create_tree (NULL, 0, 0);
1320   i = 0;
1321   curc = dest;
1322   do
1323     {
1324       all_pa_info[i].BoundingBox.X1 = curc->contours->xmin;
1325       all_pa_info[i].BoundingBox.Y1 = curc->contours->ymin;
1326       all_pa_info[i].BoundingBox.X2 = curc->contours->xmax;
1327       all_pa_info[i].BoundingBox.Y2 = curc->contours->ymax;
1328       all_pa_info[i].pa = curc;
1329       r_insert_entry (tree, (const BoxType *) &all_pa_info[i], 0);
1330       i++;
1331     }
1332   while ((curc = curc->f) != dest);
1333
1334   /* loop through the holes and put them where they belong */
1335   while ((curh = *src) != NULL)
1336     {
1337       *src = curh->next;
1338
1339       container = NULL;
1340       /* build a heap of all of the polys that the hole is inside its bounding box */
1341       heap = heap_create ();
1342       r_search (tree, (BoxType *) curh, NULL, heap_it, heap);
1343       if (heap_is_empty (heap))
1344         {
1345 #ifndef NDEBUG
1346 #ifdef DEBUG
1347           poly_dump (dest);
1348 #endif
1349 #endif
1350           poly_DelContour (&curh);
1351           error (err_bad_parm);
1352         }
1353       /* Now search the heap for the container. If there was only one item
1354        * in the heap, assume it is the container without the expense of
1355        * proving it.
1356        */
1357       pa_info = (struct polyarea_info *) heap_remove_smallest (heap);
1358       if (heap_is_empty (heap))
1359         {                       /* only one possibility it must be the right one */
1360           assert (poly_ContourInContour (pa_info->pa->contours, curh));
1361           container = pa_info->pa->contours;
1362         }
1363       else
1364         {
1365           do
1366             {
1367               if (poly_ContourInContour (pa_info->pa->contours, curh))
1368                 {
1369                   container = pa_info->pa->contours;
1370                   break;
1371                 }
1372               if (heap_is_empty (heap))
1373                 break;
1374               pa_info = (struct polyarea_info *) heap_remove_smallest (heap);
1375             }
1376           while (1);
1377         }
1378       heap_destroy (&heap);
1379       if (container == NULL)
1380         {
1381           /* bad input polygons were given */
1382 #ifndef NDEBUG
1383 #ifdef DEBUG
1384           poly_dump (dest);
1385 #endif
1386 #endif
1387           curh->next = NULL;
1388           poly_DelContour (&curh);
1389           error (err_bad_parm);
1390         }
1391       else
1392         {
1393           /* Need to check if this new hole means we need to kick out any old ones for reprocessing */
1394           while (1)
1395             {
1396               struct find_inside_info info;
1397               PLINE *prev;
1398
1399               info.want_inside = curh;
1400
1401               /* Set jump return */
1402               if (setjmp (info.jb))
1403                 {
1404                   /* Returned here! */
1405                 }
1406               else
1407                 {
1408                   info.result = NULL;
1409                   /* Rtree search, calling back a routine to longjmp back with data about any hole inside the added one */
1410                   /*   Be sure not to bother jumping back to report the main contour! */
1411                   r_search (pa_info->pa->contour_tree, (BoxType *) curh, NULL,
1412                             find_inside, &info);
1413
1414                   /* Nothing found? */
1415                   break;
1416                 }
1417
1418               /* We need to find the contour before it, so we can update its next pointer */
1419               prev = container;
1420               while (prev->next != info.result)
1421                 {
1422                   prev = prev->next;
1423                 }
1424
1425               /* Remove hole from the contour */
1426               remove_contour (pa_info->pa, prev, info.result, TRUE);
1427
1428               /* Add hole as the next on the list to be processed in this very function */
1429               info.result->next = *src;
1430               *src = info.result;
1431             }
1432           /* End check for kicked out holes */
1433
1434           /* link at front of hole list */
1435           curh->next = container->next;
1436           container->next = curh;
1437           r_insert_entry (pa_info->pa->contour_tree, (BoxType *) curh, 0);
1438
1439         }
1440     }
1441   r_destroy_tree (&tree);
1442   free (all_pa_info);
1443 }                               /* InsertHoles */
1444
1445
1446 /****************************************************************/
1447 /* routines for collecting result */
1448
1449 typedef enum
1450 {
1451   FORW, BACKW
1452 } DIRECTION;
1453
1454 /* Start Rule */
1455 typedef int (*S_Rule) (VNODE *, DIRECTION *);
1456
1457 /* Jump Rule  */
1458 typedef int (*J_Rule) (char, VNODE *, DIRECTION *);
1459
1460 static int
1461 UniteS_Rule (VNODE * cur, DIRECTION * initdir)
1462 {
1463   *initdir = FORW;
1464   return (NODE_LABEL (cur) == OUTSIDE) || (NODE_LABEL (cur) == SHARED);
1465 }
1466
1467 static int
1468 IsectS_Rule (VNODE * cur, DIRECTION * initdir)
1469 {
1470   *initdir = FORW;
1471   return (NODE_LABEL (cur) == INSIDE) || (NODE_LABEL (cur) == SHARED);
1472 }
1473
1474 static int
1475 SubS_Rule (VNODE * cur, DIRECTION * initdir)
1476 {
1477   *initdir = FORW;
1478   return (NODE_LABEL (cur) == OUTSIDE) || (NODE_LABEL (cur) == SHARED2);
1479 }
1480
1481 static int
1482 XorS_Rule (VNODE * cur, DIRECTION * initdir)
1483 {
1484   if (cur->Flags.status == INSIDE)
1485     {
1486       *initdir = BACKW;
1487       return TRUE;
1488     }
1489   if (cur->Flags.status == OUTSIDE)
1490     {
1491       *initdir = FORW;
1492       return TRUE;
1493     }
1494   return FALSE;
1495 }
1496
1497 static int
1498 IsectJ_Rule (char p, VNODE * v, DIRECTION * cdir)
1499 {
1500   assert (*cdir == FORW);
1501   return (v->Flags.status == INSIDE || v->Flags.status == SHARED);
1502 }
1503
1504 static int
1505 UniteJ_Rule (char p, VNODE * v, DIRECTION * cdir)
1506 {
1507   assert (*cdir == FORW);
1508   return (v->Flags.status == OUTSIDE || v->Flags.status == SHARED);
1509 }
1510
1511 static int
1512 XorJ_Rule (char p, VNODE * v, DIRECTION * cdir)
1513 {
1514   if (v->Flags.status == INSIDE)
1515     {
1516       *cdir = BACKW;
1517       return TRUE;
1518     }
1519   if (v->Flags.status == OUTSIDE)
1520     {
1521       *cdir = FORW;
1522       return TRUE;
1523     }
1524   return FALSE;
1525 }
1526
1527 static int
1528 SubJ_Rule (char p, VNODE * v, DIRECTION * cdir)
1529 {
1530   if (p == 'B' && v->Flags.status == INSIDE)
1531     {
1532       *cdir = BACKW;
1533       return TRUE;
1534     }
1535   if (p == 'A' && v->Flags.status == OUTSIDE)
1536     {
1537       *cdir = FORW;
1538       return TRUE;
1539     }
1540   if (v->Flags.status == SHARED2)
1541     {
1542       if (p == 'A')
1543         *cdir = FORW;
1544       else
1545         *cdir = BACKW;
1546       return TRUE;
1547     }
1548   return FALSE;
1549 }
1550
1551 /* return the edge that comes next.
1552  * if the direction is BACKW, then we return the next vertex
1553  * so that prev vertex has the flags for the edge
1554  *
1555  * returns true if an edge is found, false otherwise
1556  */
1557 static int
1558 jump (VNODE ** cur, DIRECTION * cdir, J_Rule rule)
1559 {
1560   CVCList *d, *start;
1561   VNODE *e;
1562   DIRECTION newone;
1563
1564   if (!(*cur)->cvc_prev)        /* not a cross-vertex */
1565     {
1566       if (*cdir == FORW ? (*cur)->Flags.mark : (*cur)->prev->Flags.mark)
1567         return FALSE;
1568       return TRUE;
1569     }
1570 #ifdef DEBUG_JUMP
1571   DEBUGP ("jump entering node at %$mD\n", (*cur)->point[0], (*cur)->point[1]);
1572 #endif
1573   if (*cdir == FORW)
1574     d = (*cur)->cvc_prev->prev;
1575   else
1576     d = (*cur)->cvc_next->prev;
1577   start = d;
1578   do
1579     {
1580       e = d->parent;
1581       if (d->side == 'P')
1582         e = e->prev;
1583       newone = *cdir;
1584       if (!e->Flags.mark && rule (d->poly, e, &newone))
1585         {
1586           if ((d->side == 'N' && newone == FORW) ||
1587               (d->side == 'P' && newone == BACKW))
1588             {
1589 #ifdef DEBUG_JUMP
1590               if (newone == FORW)
1591                 DEBUGP ("jump leaving node at %#mD\n",
1592                         e->next->point[0], e->next->point[1]);
1593               else
1594                 DEBUGP ("jump leaving node at %#mD\n",
1595                         e->point[0], e->point[1]);
1596 #endif
1597               *cur = d->parent;
1598               *cdir = newone;
1599               return TRUE;
1600             }
1601         }
1602     }
1603   while ((d = d->prev) != start);
1604   return FALSE;
1605 }
1606
1607 static int
1608 Gather (VNODE * start, PLINE ** result, J_Rule v_rule, DIRECTION initdir)
1609 {
1610   VNODE *cur = start, *newn;
1611   DIRECTION dir = initdir;
1612 #ifdef DEBUG_GATHER
1613   DEBUGP ("gather direction = %d\n", dir);
1614 #endif
1615   assert (*result == NULL);
1616   do
1617     {
1618       /* see where to go next */
1619       if (!jump (&cur, &dir, v_rule))
1620         break;
1621       /* add edge to polygon */
1622       if (!*result)
1623         {
1624           *result = poly_NewContour (cur->point);
1625           if (*result == NULL)
1626             return err_no_memory;
1627         }
1628       else
1629         {
1630           if ((newn = poly_CreateNode (cur->point)) == NULL)
1631             return err_no_memory;
1632           poly_InclVertex ((*result)->head.prev, newn);
1633         }
1634 #ifdef DEBUG_GATHER
1635       DEBUGP ("gather vertex at %#mD\n", cur->point[0], cur->point[1]);
1636 #endif
1637       /* Now mark the edge as included.  */
1638       newn = (dir == FORW ? cur : cur->prev);
1639       newn->Flags.mark = 1;
1640       /* for SHARED edge mark both */
1641       if (newn->shared)
1642         newn->shared->Flags.mark = 1;
1643
1644       /* Advance to the next edge.  */
1645       cur = (dir == FORW ? cur->next : cur->prev);
1646     }
1647   while (1);
1648   return err_ok;
1649 }                               /* Gather */
1650
1651 static void
1652 Collect1 (jmp_buf * e, VNODE * cur, DIRECTION dir, POLYAREA ** contours,
1653           PLINE ** holes, J_Rule j_rule)
1654 {
1655   PLINE *p = NULL;              /* start making contour */
1656   int errc = err_ok;
1657   if ((errc =
1658        Gather (dir == FORW ? cur : cur->next, &p, j_rule, dir)) != err_ok)
1659     {
1660       if (p != NULL)
1661         poly_DelContour (&p);
1662       error (errc);
1663     }
1664   if (!p)
1665     return;
1666   poly_PreContour (p, TRUE);
1667   if (p->Count > 2)
1668     {
1669 #ifdef DEBUG_GATHER
1670       DEBUGP ("adding contour with %d verticies and direction %c\n",
1671               p->Count, p->Flags.orient ? 'F' : 'B');
1672 #endif
1673       PutContour (e, p, contours, holes, NULL, NULL, NULL);
1674     }
1675   else
1676     {
1677       /* some sort of computation error ? */
1678 #ifdef DEBUG_GATHER
1679       DEBUGP ("Bad contour! Not enough points!!\n");
1680 #endif
1681       poly_DelContour (&p);
1682     }
1683 }
1684
1685 static void
1686 Collect (jmp_buf * e, PLINE * a, POLYAREA ** contours, PLINE ** holes,
1687          S_Rule s_rule, J_Rule j_rule)
1688 {
1689   VNODE *cur, *other;
1690   DIRECTION dir;
1691
1692   cur = &a->head;
1693   do
1694     {
1695       if (s_rule (cur, &dir) && cur->Flags.mark == 0)
1696         Collect1 (e, cur, dir, contours, holes, j_rule);
1697       other = cur;
1698       if ((other->cvc_prev && jump (&other, &dir, j_rule)))
1699         Collect1 (e, other, dir, contours, holes, j_rule);
1700     }
1701   while ((cur = cur->next) != &a->head);
1702 }                               /* Collect */
1703
1704
1705 static int
1706 cntr_Collect (jmp_buf * e, PLINE ** A, POLYAREA ** contours, PLINE ** holes,
1707               int action, POLYAREA * owner, POLYAREA * parent,
1708               PLINE * parent_contour)
1709 {
1710   PLINE *tmprev;
1711
1712   if ((*A)->Flags.status == ISECTED)
1713     {
1714       switch (action)
1715         {
1716         case PBO_UNITE:
1717           Collect (e, *A, contours, holes, UniteS_Rule, UniteJ_Rule);
1718           break;
1719         case PBO_ISECT:
1720           Collect (e, *A, contours, holes, IsectS_Rule, IsectJ_Rule);
1721           break;
1722         case PBO_XOR:
1723           Collect (e, *A, contours, holes, XorS_Rule, XorJ_Rule);
1724           break;
1725         case PBO_SUB:
1726           Collect (e, *A, contours, holes, SubS_Rule, SubJ_Rule);
1727           break;
1728         };
1729     }
1730   else
1731     {
1732       switch (action)
1733         {
1734         case PBO_ISECT:
1735           if ((*A)->Flags.status == INSIDE)
1736             {
1737               tmprev = *A;
1738               /* disappear this contour (rtree entry removed in PutContour) */
1739               *A = tmprev->next;
1740               tmprev->next = NULL;
1741               PutContour (e, tmprev, contours, holes, owner, NULL, NULL);
1742               return TRUE;
1743             }
1744           break;
1745         case PBO_XOR:
1746           if ((*A)->Flags.status == INSIDE)
1747             {
1748               tmprev = *A;
1749               /* disappear this contour (rtree entry removed in PutContour) */
1750               *A = tmprev->next;
1751               tmprev->next = NULL;
1752               poly_InvContour (tmprev);
1753               PutContour (e, tmprev, contours, holes, owner, NULL, NULL);
1754               return TRUE;
1755             }
1756           /* break; *//* Fall through (I think this is correct! pcjc2) */
1757         case PBO_UNITE:
1758         case PBO_SUB:
1759           if ((*A)->Flags.status == OUTSIDE)
1760             {
1761               tmprev = *A;
1762               /* disappear this contour (rtree entry removed in PutContour) */
1763               *A = tmprev->next;
1764               tmprev->next = NULL;
1765               PutContour (e, tmprev, contours, holes, owner, parent,
1766                           parent_contour);
1767               return TRUE;
1768             }
1769           break;
1770         }
1771     }
1772   return FALSE;
1773 }                               /* cntr_Collect */
1774
1775 static void
1776 M_B_AREA_Collect (jmp_buf * e, POLYAREA * bfst, POLYAREA ** contours,
1777                   PLINE ** holes, int action)
1778 {
1779   POLYAREA *b = bfst;
1780   PLINE **cur, **next, *tmp;
1781
1782   assert (b != NULL);
1783   do
1784     {
1785       for (cur = &b->contours; *cur != NULL; cur = next)
1786         {
1787           next = &((*cur)->next);
1788           if ((*cur)->Flags.status == ISECTED)
1789             continue;
1790
1791           if ((*cur)->Flags.status == INSIDE)
1792             switch (action)
1793               {
1794               case PBO_XOR:
1795               case PBO_SUB:
1796                 poly_InvContour (*cur);
1797               case PBO_ISECT:
1798                 tmp = *cur;
1799                 *cur = tmp->next;
1800                 next = cur;
1801                 tmp->next = NULL;
1802                 tmp->Flags.status = UNKNWN;
1803                 PutContour (e, tmp, contours, holes, b, NULL, NULL);
1804                 break;
1805               case PBO_UNITE:
1806                 break;          /* nothing to do - already included */
1807               }
1808           else if ((*cur)->Flags.status == OUTSIDE)
1809             switch (action)
1810               {
1811               case PBO_XOR:
1812               case PBO_UNITE:
1813                 /* include */
1814                 tmp = *cur;
1815                 *cur = tmp->next;
1816                 next = cur;
1817                 tmp->next = NULL;
1818                 tmp->Flags.status = UNKNWN;
1819                 PutContour (e, tmp, contours, holes, b, NULL, NULL);
1820                 break;
1821               case PBO_ISECT:
1822               case PBO_SUB:
1823                 break;          /* do nothing, not included */
1824               }
1825         }
1826     }
1827   while ((b = b->f) != bfst);
1828 }
1829
1830
1831 static inline int
1832 contour_is_first (POLYAREA * a, PLINE * cur)
1833 {
1834   return (a->contours == cur);
1835 }
1836
1837
1838 static inline int
1839 contour_is_last (PLINE * cur)
1840 {
1841   return (cur->next == NULL);
1842 }
1843
1844
1845 static inline void
1846 remove_polyarea (POLYAREA ** list, POLYAREA * piece)
1847 {
1848   /* If this item was the start of the list, advance that pointer */
1849   if (*list == piece)
1850     *list = (*list)->f;
1851
1852   /* But reset it to NULL if it wraps around and hits us again */
1853   if (*list == piece)
1854     *list = NULL;
1855
1856   piece->b->f = piece->f;
1857   piece->f->b = piece->b;
1858   piece->f = piece->b = piece;
1859 }
1860
1861 static void
1862 M_POLYAREA_separate_isected (jmp_buf * e, POLYAREA ** pieces,
1863                              PLINE ** holes, PLINE ** isected)
1864 {
1865   POLYAREA *a = *pieces;
1866   POLYAREA *anext;
1867   PLINE *curc, *next, *prev;
1868   int finished;
1869
1870   if (a == NULL)
1871     return;
1872
1873   /* TODO: STASH ENOUGH INFORMATION EARLIER ON, SO WE CAN REMOVE THE INTERSECTED
1874      CONTOURS WITHOUT HAVING TO WALK THE FULL DATA-STRUCTURE LOOKING FOR THEM. */
1875
1876   do
1877     {
1878       int hole_contour = 0;
1879       int is_outline = 1;
1880
1881       anext = a->f;
1882       finished = (anext == *pieces);
1883
1884       prev = NULL;
1885       for (curc = a->contours; curc != NULL; curc = next, is_outline = 0)
1886         {
1887           int is_first = contour_is_first (a, curc);
1888           int is_last = contour_is_last (curc);
1889           int isect_contour = (curc->Flags.status == ISECTED);
1890
1891           next = curc->next;
1892
1893           if (isect_contour || hole_contour)
1894             {
1895
1896               /* Reset the intersection flags, since we keep these pieces */
1897               if (curc->Flags.status != ISECTED)
1898                 curc->Flags.status = UNKNWN;
1899
1900               remove_contour (a, prev, curc, !(is_first && is_last));
1901
1902               if (isect_contour)
1903                 {
1904                   /* Link into the list of intersected contours */
1905                   curc->next = *isected;
1906                   *isected = curc;
1907                 }
1908               else if (hole_contour)
1909                 {
1910                   /* Link into the list of holes */
1911                   curc->next = *holes;
1912                   *holes = curc;
1913                 }
1914               else
1915                 {
1916                   assert (0);
1917                 }
1918
1919               if (is_first && is_last)
1920                 {
1921                   remove_polyarea (pieces, a);
1922                   poly_Free (&a);       /* NB: Sets a to NULL */
1923                 }
1924
1925             }
1926           else
1927             {
1928               /* Note the item we just didn't delete as the next
1929                  candidate for having its "next" pointer adjusted.
1930                  Saves walking the contour list when we delete one. */
1931               prev = curc;
1932             }
1933
1934           /* If we move or delete an outer contour, we need to move any holes
1935              we wish to keep within that contour to the holes list. */
1936           if (is_outline && isect_contour)
1937             hole_contour = 1;
1938
1939         }
1940
1941       /* If we deleted all the pieces of the polyarea, *pieces is NULL */
1942     }
1943   while ((a = anext), *pieces != NULL && !finished);
1944 }
1945
1946
1947 struct find_inside_m_pa_info
1948 {
1949   jmp_buf jb;
1950   POLYAREA *want_inside;
1951   PLINE *result;
1952 };
1953
1954 static int
1955 find_inside_m_pa (const BoxType * b, void *cl)
1956 {
1957   struct find_inside_m_pa_info *info = (struct find_inside_m_pa_info *) cl;
1958   PLINE *check = (PLINE *) b;
1959   /* Don't report for the main contour */
1960   if (check->Flags.orient == PLF_DIR)
1961     return 0;
1962   /* Don't look at contours marked as being intersected */
1963   if (check->Flags.status == ISECTED)
1964     return 0;
1965   if (cntr_in_M_POLYAREA (check, info->want_inside, FALSE))
1966     {
1967       info->result = check;
1968       longjmp (info->jb, 1);
1969     }
1970   return 0;
1971 }
1972
1973
1974 static void
1975 M_POLYAREA_update_primary (jmp_buf * e, POLYAREA ** pieces,
1976                            PLINE ** holes, int action, POLYAREA * bpa)
1977 {
1978   POLYAREA *a = *pieces;
1979   POLYAREA *b;
1980   POLYAREA *anext;
1981   PLINE *curc, *next, *prev;
1982   BoxType box;
1983   /* int inv_inside = 0; */
1984   int del_inside = 0;
1985   int del_outside = 0;
1986   int finished;
1987
1988   if (a == NULL)
1989     return;
1990
1991   switch (action)
1992     {
1993     case PBO_ISECT:
1994       del_outside = 1;
1995       break;
1996     case PBO_UNITE:
1997     case PBO_SUB:
1998       del_inside = 1;
1999       break;
2000     case PBO_XOR:               /* NOT IMPLEMENTED OR USED */
2001       /* inv_inside = 1; */
2002       assert (0);
2003       break;
2004     }
2005
2006   box = *((BoxType *) bpa->contours);
2007   b = bpa;
2008   while ((b = b->f) != bpa)
2009     {
2010       BoxType *b_box = (BoxType *) b->contours;
2011       MAKEMIN (box.X1, b_box->X1);
2012       MAKEMIN (box.Y1, b_box->Y1);
2013       MAKEMAX (box.X2, b_box->X2);
2014       MAKEMAX (box.Y2, b_box->Y2);
2015     }
2016
2017   if (del_inside)
2018     {
2019
2020       do
2021         {
2022           anext = a->f;
2023           finished = (anext == *pieces);
2024
2025           /* Test the outer contour first, as we may need to remove all children */
2026
2027           /* We've not yet split intersected contours out, just ignore them */
2028           if (a->contours->Flags.status != ISECTED &&
2029               /* Pre-filter on bounding box */
2030               ((a->contours->xmin >= box.X1) && (a->contours->ymin >= box.Y1)
2031                && (a->contours->xmax <= box.X2)
2032                && (a->contours->ymax <= box.Y2)) &&
2033               /* Then test properly */
2034               cntr_in_M_POLYAREA (a->contours, bpa, FALSE))
2035             {
2036
2037               /* Delete this contour, all children -> holes queue */
2038
2039               /* Delete the outer contour */
2040               curc = a->contours;
2041               remove_contour (a, NULL, curc, FALSE);    /* Rtree deleted in poly_Free below */
2042               /* a->contours now points to the remaining holes */
2043               poly_DelContour (&curc);
2044
2045               if (a->contours != NULL)
2046                 {
2047                   /* Find the end of the list of holes */
2048                   curc = a->contours;
2049                   while (curc->next != NULL)
2050                     curc = curc->next;
2051
2052                   /* Take the holes and prepend to the holes queue */
2053                   curc->next = *holes;
2054                   *holes = a->contours;
2055                   a->contours = NULL;
2056                 }
2057
2058               remove_polyarea (pieces, a);
2059               poly_Free (&a);   /* NB: Sets a to NULL */
2060
2061               continue;
2062             }
2063
2064           /* Loop whilst we find INSIDE contours to delete */
2065           while (1)
2066             {
2067               struct find_inside_m_pa_info info;
2068               PLINE *prev;
2069
2070               info.want_inside = bpa;
2071
2072               /* Set jump return */
2073               if (setjmp (info.jb))
2074                 {
2075                   /* Returned here! */
2076                 }
2077               else
2078                 {
2079                   info.result = NULL;
2080                   /* r-tree search, calling back a routine to longjmp back with
2081                    * data about any hole inside the B polygon.
2082                    * NB: Does not jump back to report the main contour!
2083                    */
2084                   r_search (a->contour_tree, &box, NULL, find_inside_m_pa,
2085                             &info);
2086
2087                   /* Nothing found? */
2088                   break;
2089                 }
2090
2091               /* We need to find the contour before it, so we can update its next pointer */
2092               prev = a->contours;
2093               while (prev->next != info.result)
2094                 {
2095                   prev = prev->next;
2096                 }
2097
2098               /* Remove hole from the contour */
2099               remove_contour (a, prev, info.result, TRUE);
2100               poly_DelContour (&info.result);
2101             }
2102           /* End check for deleted holes */
2103
2104           /* If we deleted all the pieces of the polyarea, *pieces is NULL */
2105         }
2106       while ((a = anext), *pieces != NULL && !finished);
2107
2108       return;
2109     }
2110   else
2111     {
2112       /* This path isn't optimised for speed */
2113     }
2114
2115   do
2116     {
2117       int hole_contour = 0;
2118       int is_outline = 1;
2119
2120       anext = a->f;
2121       finished = (anext == *pieces);
2122
2123       prev = NULL;
2124       for (curc = a->contours; curc != NULL; curc = next, is_outline = 0)
2125         {
2126           int is_first = contour_is_first (a, curc);
2127           int is_last = contour_is_last (curc);
2128           int del_contour = 0;
2129
2130           next = curc->next;
2131
2132           if (del_outside)
2133             del_contour = curc->Flags.status != ISECTED &&
2134               !cntr_in_M_POLYAREA (curc, bpa, FALSE);
2135
2136           /* Skip intersected contours */
2137           if (curc->Flags.status == ISECTED)
2138             {
2139               prev = curc;
2140               continue;
2141             }
2142
2143           /* Reset the intersection flags, since we keep these pieces */
2144           curc->Flags.status = UNKNWN;
2145
2146           if (del_contour || hole_contour)
2147             {
2148
2149               remove_contour (a, prev, curc, !(is_first && is_last));
2150
2151               if (del_contour)
2152                 {
2153                   /* Delete the contour */
2154                   poly_DelContour (&curc);      /* NB: Sets curc to NULL */
2155                 }
2156               else if (hole_contour)
2157                 {
2158                   /* Link into the list of holes */
2159                   curc->next = *holes;
2160                   *holes = curc;
2161                 }
2162               else
2163                 {
2164                   assert (0);
2165                 }
2166
2167               if (is_first && is_last)
2168                 {
2169                   remove_polyarea (pieces, a);
2170                   poly_Free (&a);       /* NB: Sets a to NULL */
2171                 }
2172
2173             }
2174           else
2175             {
2176               /* Note the item we just didn't delete as the next
2177                  candidate for having its "next" pointer adjusted.
2178                  Saves walking the contour list when we delete one. */
2179               prev = curc;
2180             }
2181
2182           /* If we move or delete an outer contour, we need to move any holes
2183              we wish to keep within that contour to the holes list. */
2184           if (is_outline && del_contour)
2185             hole_contour = 1;
2186
2187         }
2188
2189       /* If we deleted all the pieces of the polyarea, *pieces is NULL */
2190     }
2191   while ((a = anext), *pieces != NULL && !finished);
2192 }
2193
2194 static void
2195 M_POLYAREA_Collect_separated (jmp_buf * e, PLINE * afst, POLYAREA ** contours,
2196                               PLINE ** holes, int action, BOOLp maybe)
2197 {
2198   PLINE **cur, **next;
2199
2200   for (cur = &afst; *cur != NULL; cur = next)
2201     {
2202       next = &((*cur)->next);
2203       /* if we disappear a contour, don't advance twice */
2204       if (cntr_Collect (e, cur, contours, holes, action, NULL, NULL, NULL))
2205         next = cur;
2206     }
2207 }
2208
2209 static void
2210 M_POLYAREA_Collect (jmp_buf * e, POLYAREA * afst, POLYAREA ** contours,
2211                     PLINE ** holes, int action, BOOLp maybe)
2212 {
2213   POLYAREA *a = afst;
2214   POLYAREA *parent = NULL;      /* Quiet compiler warning */
2215   PLINE **cur, **next, *parent_contour;
2216
2217   assert (a != NULL);
2218   while ((a = a->f) != afst);
2219   /* now the non-intersect parts are collected in temp/holes */
2220   do
2221     {
2222       if (maybe && a->contours->Flags.status != ISECTED)
2223         parent_contour = a->contours;
2224       else
2225         parent_contour = NULL;
2226
2227       /* Take care of the first contour - so we know if we
2228        * can shortcut reparenting some of its children
2229        */
2230       cur = &a->contours;
2231       if (*cur != NULL)
2232         {
2233           next = &((*cur)->next);
2234           /* if we disappear a contour, don't advance twice */
2235           if (cntr_Collect (e, cur, contours, holes, action, a, NULL, NULL))
2236             {
2237               parent = (*contours)->b;  /* InsCntr inserts behind the head */
2238               next = cur;
2239             }
2240           else
2241             parent = a;
2242           cur = next;
2243         }
2244       for (; *cur != NULL; cur = next)
2245         {
2246           next = &((*cur)->next);
2247           /* if we disappear a contour, don't advance twice */
2248           if (cntr_Collect (e, cur, contours, holes, action, a, parent,
2249                             parent_contour))
2250             next = cur;
2251         }
2252     }
2253   while ((a = a->f) != afst);
2254 }
2255
2256 /* determine if two polygons touch or overlap */
2257 BOOLp
2258 Touching (POLYAREA * a, POLYAREA * b)
2259 {
2260   jmp_buf e;
2261   int code;
2262
2263   if ((code = setjmp (e)) == 0)
2264     {
2265 #ifdef DEBUG
2266       if (!poly_Valid (a))
2267         return -1;
2268       if (!poly_Valid (b))
2269         return -1;
2270 #endif
2271       M_POLYAREA_intersect (&e, a, b, false);
2272
2273       if (M_POLYAREA_label (a, b, TRUE))
2274         return TRUE;
2275       if (M_POLYAREA_label (b, a, TRUE))
2276         return TRUE;
2277     }
2278   else if (code == TOUCHES)
2279     return TRUE;
2280   return FALSE;
2281 }
2282
2283 /* the main clipping routines */
2284 int
2285 poly_Boolean (const POLYAREA * a_org, const POLYAREA * b_org,
2286               POLYAREA ** res, int action)
2287 {
2288   POLYAREA *a = NULL, *b = NULL;
2289
2290   if (!poly_M_Copy0 (&a, a_org) || !poly_M_Copy0 (&b, b_org))
2291     return err_no_memory;
2292
2293   return poly_Boolean_free (a, b, res, action);
2294 }                               /* poly_Boolean */
2295
2296 /* just like poly_Boolean but frees the input polys */
2297 int
2298 poly_Boolean_free (POLYAREA * ai, POLYAREA * bi, POLYAREA ** res, int action)
2299 {
2300   POLYAREA *a = ai, *b = bi;
2301   PLINE *a_isected = NULL;
2302   PLINE *p, *holes = NULL;
2303   jmp_buf e;
2304   int code;
2305
2306   *res = NULL;
2307
2308   if (!a)
2309     {
2310       switch (action)
2311         {
2312         case PBO_XOR:
2313         case PBO_UNITE:
2314           *res = bi;
2315           return err_ok;
2316         case PBO_SUB:
2317         case PBO_ISECT:
2318           if (b != NULL)
2319             poly_Free (&b);
2320           return err_ok;
2321         }
2322     }
2323   if (!b)
2324     {
2325       switch (action)
2326         {
2327         case PBO_SUB:
2328         case PBO_XOR:
2329         case PBO_UNITE:
2330           *res = ai;
2331           return err_ok;
2332         case PBO_ISECT:
2333           if (a != NULL)
2334             poly_Free (&a);
2335           return err_ok;
2336         }
2337     }
2338
2339   if ((code = setjmp (e)) == 0)
2340     {
2341 #ifdef DEBUG
2342       assert (poly_Valid (a));
2343       assert (poly_Valid (b));
2344 #endif
2345
2346       /* intersect needs to make a list of the contours in a and b which are intersected */
2347       M_POLYAREA_intersect (&e, a, b, TRUE);
2348
2349       /* We could speed things up a lot here if we only processed the relevant contours */
2350       /* NB: Relevant parts of a are labeled below */
2351       M_POLYAREA_label (b, a, FALSE);
2352
2353       *res = a;
2354       M_POLYAREA_update_primary (&e, res, &holes, action, b);
2355       M_POLYAREA_separate_isected (&e, res, &holes, &a_isected);
2356       M_POLYAREA_label_separated (a_isected, b, FALSE);
2357       M_POLYAREA_Collect_separated (&e, a_isected, res, &holes, action,
2358                                     FALSE);
2359       M_B_AREA_Collect (&e, b, res, &holes, action);
2360       poly_Free (&b);
2361
2362       /* free a_isected */
2363       while ((p = a_isected) != NULL)
2364         {
2365           a_isected = p->next;
2366           poly_DelContour (&p);
2367         }
2368
2369       InsertHoles (&e, *res, &holes);
2370     }
2371   /* delete holes if any left */
2372   while ((p = holes) != NULL)
2373     {
2374       holes = p->next;
2375       poly_DelContour (&p);
2376     }
2377
2378   if (code)
2379     {
2380       poly_Free (res);
2381       return code;
2382     }
2383   assert (!*res || poly_Valid (*res));
2384   return code;
2385 }                               /* poly_Boolean_free */
2386
2387 static void
2388 clear_marks (POLYAREA * p)
2389 {
2390   POLYAREA *n = p;
2391   PLINE *c;
2392   VNODE *v;
2393
2394   do
2395     {
2396       for (c = n->contours; c; c = c->next)
2397         {
2398           v = &c->head;
2399           do
2400             {
2401               v->Flags.mark = 0;
2402             }
2403           while ((v = v->next) != &c->head);
2404         }
2405     }
2406   while ((n = n->f) != p);
2407 }
2408
2409 /* compute the intersection and subtraction (divides "a" into two pieces)
2410  * and frees the input polys. This assumes that bi is a single simple polygon.
2411  */
2412 int
2413 poly_AndSubtract_free (POLYAREA * ai, POLYAREA * bi,
2414                        POLYAREA ** aandb, POLYAREA ** aminusb)
2415 {
2416   POLYAREA *a = ai, *b = bi;
2417   PLINE *p, *holes = NULL;
2418   jmp_buf e;
2419   int code;
2420
2421   *aandb = NULL;
2422   *aminusb = NULL;
2423
2424   if ((code = setjmp (e)) == 0)
2425     {
2426
2427 #ifdef DEBUG
2428       if (!poly_Valid (a))
2429         return -1;
2430       if (!poly_Valid (b))
2431         return -1;
2432 #endif
2433       M_POLYAREA_intersect (&e, a, b, TRUE);
2434
2435       M_POLYAREA_label (a, b, FALSE);
2436       M_POLYAREA_label (b, a, FALSE);
2437
2438       M_POLYAREA_Collect (&e, a, aandb, &holes, PBO_ISECT, FALSE);
2439       InsertHoles (&e, *aandb, &holes);
2440       assert (poly_Valid (*aandb));
2441       /* delete holes if any left */
2442       while ((p = holes) != NULL)
2443         {
2444           holes = p->next;
2445           poly_DelContour (&p);
2446         }
2447       holes = NULL;
2448       clear_marks (a);
2449       clear_marks (b);
2450       M_POLYAREA_Collect (&e, a, aminusb, &holes, PBO_SUB, FALSE);
2451       InsertHoles (&e, *aminusb, &holes);
2452       poly_Free (&a);
2453       poly_Free (&b);
2454       assert (poly_Valid (*aminusb));
2455     }
2456   /* delete holes if any left */
2457   while ((p = holes) != NULL)
2458     {
2459       holes = p->next;
2460       poly_DelContour (&p);
2461     }
2462
2463
2464   if (code)
2465     {
2466       poly_Free (aandb);
2467       poly_Free (aminusb);
2468       return code;
2469     }
2470   assert (!*aandb || poly_Valid (*aandb));
2471   assert (!*aminusb || poly_Valid (*aminusb));
2472   return code;
2473 }                               /* poly_AndSubtract_free */
2474
2475 static inline int
2476 cntrbox_pointin (PLINE * c, Vector p)
2477 {
2478   return (p[0] >= c->xmin && p[1] >= c->ymin &&
2479           p[0] <= c->xmax && p[1] <= c->ymax);
2480
2481 }
2482
2483 static inline int
2484 node_neighbours (VNODE * a, VNODE * b)
2485 {
2486   return (a == b) || (a->next == b) || (b->next == a) || (a->next == b->next);
2487 }
2488
2489 VNODE *
2490 poly_CreateNode (Vector v)
2491 {
2492   VNODE *res;
2493   Coord *c;
2494
2495   assert (v);
2496   res = (VNODE *) calloc (1, sizeof (VNODE));
2497   if (res == NULL)
2498     return NULL;
2499   // bzero (res, sizeof (VNODE) - sizeof(Vector));
2500   c = res->point;
2501   *c++ = *v++;
2502   *c = *v;
2503   return res;
2504 }
2505
2506 void
2507 poly_IniContour (PLINE * c)
2508 {
2509   if (c == NULL)
2510     return;
2511   /* bzero (c, sizeof(PLINE)); */
2512   c->head.next = c->head.prev = &c->head;
2513   c->xmin = c->ymin = COORD_MAX;
2514   c->xmax = c->ymax = -COORD_MAX - 1;
2515   c->is_round = FALSE;
2516   c->cx = 0;
2517   c->cy = 0;
2518   c->radius = 0;
2519 }
2520
2521 PLINE *
2522 poly_NewContour (Vector v)
2523 {
2524   PLINE *res;
2525
2526   res = (PLINE *) calloc (1, sizeof (PLINE));
2527   if (res == NULL)
2528     return NULL;
2529
2530   poly_IniContour (res);
2531
2532   if (v != NULL)
2533     {
2534       Vcopy (res->head.point, v);
2535       cntrbox_adjust (res, v);
2536     }
2537
2538   return res;
2539 }
2540
2541 void
2542 poly_ClrContour (PLINE * c)
2543 {
2544   VNODE *cur;
2545
2546   assert (c != NULL);
2547   while ((cur = c->head.next) != &c->head)
2548     {
2549       poly_ExclVertex (cur);
2550       free (cur);
2551     }
2552   poly_IniContour (c);
2553 }
2554
2555 void
2556 poly_DelContour (PLINE ** c)
2557 {
2558   VNODE *cur, *prev;
2559
2560   if (*c == NULL)
2561     return;
2562   for (cur = (*c)->head.prev; cur != &(*c)->head; cur = prev)
2563     {
2564       prev = cur->prev;
2565       if (cur->cvc_next != NULL)
2566         {
2567           free (cur->cvc_next);
2568           free (cur->cvc_prev);
2569         }
2570       free (cur);
2571     }
2572   if ((*c)->head.cvc_next != NULL)
2573     {
2574       free ((*c)->head.cvc_next);
2575       free ((*c)->head.cvc_prev);
2576     }
2577   /* FIXME -- strict aliasing violation.  */
2578   if ((*c)->tree)
2579     {
2580       rtree_t *r = (*c)->tree;
2581       r_destroy_tree (&r);
2582     }
2583   free (*c), *c = NULL;
2584 }
2585
2586 void
2587 poly_PreContour (PLINE * C, BOOLp optimize)
2588 {
2589   double area = 0;
2590   VNODE *p, *c;
2591   Vector p1, p2;
2592
2593   assert (C != NULL);
2594
2595   if (optimize)
2596     {
2597       for (c = (p = &C->head)->next; c != &C->head; c = (p = c)->next)
2598         {
2599           /* if the previous node is on the same line with this one, we should remove it */
2600           Vsub2 (p1, c->point, p->point);
2601           Vsub2 (p2, c->next->point, c->point);
2602           /* If the product below is zero then
2603            * the points on either side of c
2604            * are on the same line!
2605            * So, remove the point c
2606            */
2607
2608           if (vect_det2 (p1, p2) == 0)
2609             {
2610               poly_ExclVertex (c);
2611               free (c);
2612               c = p;
2613             }
2614         }
2615     }
2616   C->Count = 0;
2617   C->xmin = C->xmax = C->head.point[0];
2618   C->ymin = C->ymax = C->head.point[1];
2619
2620   p = (c = &C->head)->prev;
2621   if (c != p)
2622     {
2623       do
2624         {
2625           /* calculate area for orientation */
2626           area +=
2627             (double) (p->point[0] - c->point[0]) * (p->point[1] +
2628                                                     c->point[1]);
2629           cntrbox_adjust (C, c->point);
2630           C->Count++;
2631         }
2632       while ((c = (p = c)->next) != &C->head);
2633     }
2634   C->area = ABS (area);
2635   if (C->Count > 2)
2636     C->Flags.orient = ((area < 0) ? PLF_INV : PLF_DIR);
2637   C->tree = (rtree_t *)make_edge_tree (C);
2638 }                               /* poly_PreContour */
2639
2640 static int
2641 flip_cb (const BoxType * b, void *cl)
2642 {
2643   struct seg *s = (struct seg *) b;
2644   s->v = s->v->prev;
2645   return 1;
2646 }
2647
2648 void
2649 poly_InvContour (PLINE * c)
2650 {
2651   VNODE *cur, *next;
2652 #ifndef NDEBUG
2653   int r;
2654 #endif
2655
2656   assert (c != NULL);
2657   cur = &c->head;
2658   do
2659     {
2660       next = cur->next;
2661       cur->next = cur->prev;
2662       cur->prev = next;
2663       /* fix the segment tree */
2664     }
2665   while ((cur = next) != &c->head);
2666   c->Flags.orient ^= 1;
2667   if (c->tree)
2668     {
2669 #ifndef NDEBUG
2670       r =
2671 #endif
2672         r_search (c->tree, NULL, NULL, flip_cb, NULL);
2673 #ifndef NDEBUG
2674       assert (r == c->Count);
2675 #endif
2676     }
2677 }
2678
2679 void
2680 poly_ExclVertex (VNODE * node)
2681 {
2682   assert (node != NULL);
2683   if (node->cvc_next)
2684     {
2685       free (node->cvc_next);
2686       free (node->cvc_prev);
2687     }
2688   node->prev->next = node->next;
2689   node->next->prev = node->prev;
2690 }
2691
2692 void
2693 poly_InclVertex (VNODE * after, VNODE * node)
2694 {
2695   double a, b;
2696   assert (after != NULL);
2697   assert (node != NULL);
2698
2699   node->prev = after;
2700   node->next = after->next;
2701   after->next = after->next->prev = node;
2702   /* remove points on same line */
2703   if (node->prev->prev == node)
2704     return;                     /* we don't have 3 points in the poly yet */
2705   a = (node->point[1] - node->prev->prev->point[1]);
2706   a *= (node->prev->point[0] - node->prev->prev->point[0]);
2707   b = (node->point[0] - node->prev->prev->point[0]);
2708   b *= (node->prev->point[1] - node->prev->prev->point[1]);
2709   if (fabs (a - b) < EPSILON)
2710     {
2711       VNODE *t = node->prev;
2712       t->prev->next = node;
2713       node->prev = t->prev;
2714       free (t);
2715     }
2716 }
2717
2718 BOOLp
2719 poly_CopyContour (PLINE ** dst, PLINE * src)
2720 {
2721   VNODE *cur, *newnode;
2722
2723   assert (src != NULL);
2724   *dst = NULL;
2725   *dst = poly_NewContour (src->head.point);
2726   if (*dst == NULL)
2727     return FALSE;
2728
2729   (*dst)->Count = src->Count;
2730   (*dst)->Flags.orient = src->Flags.orient;
2731   (*dst)->xmin = src->xmin, (*dst)->xmax = src->xmax;
2732   (*dst)->ymin = src->ymin, (*dst)->ymax = src->ymax;
2733   (*dst)->area = src->area;
2734
2735   for (cur = src->head.next; cur != &src->head; cur = cur->next)
2736     {
2737       if ((newnode = poly_CreateNode (cur->point)) == NULL)
2738         return FALSE;
2739       // newnode->Flags = cur->Flags;
2740       poly_InclVertex ((*dst)->head.prev, newnode);
2741     }
2742   (*dst)->tree = (rtree_t *)make_edge_tree (*dst);
2743   return TRUE;
2744 }
2745
2746 /**********************************************************************/
2747 /* polygon routines */
2748
2749 BOOLp
2750 poly_Copy0 (POLYAREA ** dst, const POLYAREA * src)
2751 {
2752   *dst = NULL;
2753   if (src != NULL)
2754     *dst = (POLYAREA *)calloc (1, sizeof (POLYAREA));
2755   if (*dst == NULL)
2756     return FALSE;
2757   (*dst)->contour_tree = r_create_tree (NULL, 0, 0);
2758
2759   return poly_Copy1 (*dst, src);
2760 }
2761
2762 BOOLp
2763 poly_Copy1 (POLYAREA * dst, const POLYAREA * src)
2764 {
2765   PLINE *cur, **last = &dst->contours;
2766
2767   *last = NULL;
2768   dst->f = dst->b = dst;
2769
2770   for (cur = src->contours; cur != NULL; cur = cur->next)
2771     {
2772       if (!poly_CopyContour (last, cur))
2773         return FALSE;
2774       r_insert_entry (dst->contour_tree, (BoxType *) * last, 0);
2775       last = &(*last)->next;
2776     }
2777   return TRUE;
2778 }
2779
2780 void
2781 poly_M_Incl (POLYAREA ** list, POLYAREA * a)
2782 {
2783   if (*list == NULL)
2784     a->f = a->b = a, *list = a;
2785   else
2786     {
2787       a->f = *list;
2788       a->b = (*list)->b;
2789       (*list)->b = (*list)->b->f = a;
2790     }
2791 }
2792
2793 BOOLp
2794 poly_M_Copy0 (POLYAREA ** dst, const POLYAREA * srcfst)
2795 {
2796   const POLYAREA *src = srcfst;
2797   POLYAREA *di;
2798
2799   *dst = NULL;
2800   if (src == NULL)
2801     return FALSE;
2802   do
2803     {
2804       if ((di = poly_Create ()) == NULL || !poly_Copy1 (di, src))
2805         return FALSE;
2806       poly_M_Incl (dst, di);
2807     }
2808   while ((src = src->f) != srcfst);
2809   return TRUE;
2810 }
2811
2812 BOOLp
2813 poly_InclContour (POLYAREA * p, PLINE * c)
2814 {
2815   PLINE *tmp;
2816
2817   if ((c == NULL) || (p == NULL))
2818     return FALSE;
2819   if (c->Flags.orient == PLF_DIR)
2820     {
2821       if (p->contours != NULL)
2822         return FALSE;
2823       p->contours = c;
2824     }
2825   else
2826     {
2827       if (p->contours == NULL)
2828         return FALSE;
2829       /* link at front of hole list */
2830       tmp = p->contours->next;
2831       p->contours->next = c;
2832       c->next = tmp;
2833     }
2834   r_insert_entry (p->contour_tree, (BoxType *) c, 0);
2835   return TRUE;
2836 }
2837
2838 typedef struct pip
2839 {
2840   int f;
2841   Vector p;
2842   jmp_buf env;
2843 } pip;
2844
2845
2846 static int
2847 crossing (const BoxType * b, void *cl)
2848 {
2849   struct seg *s = (struct seg *) b;
2850   struct pip *p = (struct pip *) cl;
2851
2852   if (s->v->point[1] <= p->p[1])
2853     {
2854       if (s->v->next->point[1] > p->p[1])
2855         {
2856           Vector v1, v2;
2857           long long cross;
2858           Vsub2 (v1, s->v->next->point, s->v->point);
2859           Vsub2 (v2, p->p, s->v->point);
2860           cross = (long long) v1[0] * v2[1] - (long long) v2[0] * v1[1];
2861           if (cross == 0)
2862             {
2863               p->f = 1;
2864               longjmp (p->env, 1);
2865             }
2866           if (cross > 0)
2867             p->f += 1;
2868         }
2869     }
2870   else
2871     {
2872       if (s->v->next->point[1] <= p->p[1])
2873         {
2874           Vector v1, v2;
2875           long long cross;
2876           Vsub2 (v1, s->v->next->point, s->v->point);
2877           Vsub2 (v2, p->p, s->v->point);
2878           cross = (long long) v1[0] * v2[1] - (long long) v2[0] * v1[1];
2879           if (cross == 0)
2880             {
2881               p->f = 1;
2882               longjmp (p->env, 1);
2883             }
2884           if (cross < 0)
2885             p->f -= 1;
2886         }
2887     }
2888   return 1;
2889 }
2890
2891 int
2892 poly_InsideContour (PLINE * c, Vector p)
2893 {
2894   struct pip info;
2895   BoxType ray;
2896
2897   if (!cntrbox_pointin (c, p))
2898     return FALSE;
2899   info.f = 0;
2900   info.p[0] = ray.X1 = p[0];
2901   info.p[1] = ray.Y1 = p[1];
2902   ray.X2 = COORD_MAX;
2903   ray.Y2 = p[1] + 1;
2904   if (setjmp (info.env) == 0)
2905     r_search (c->tree, &ray, NULL, crossing, &info);
2906   return info.f;
2907 }
2908
2909 BOOLp
2910 poly_CheckInside (POLYAREA * p, Vector v0)
2911 {
2912   PLINE *cur;
2913
2914   if ((p == NULL) || (v0 == NULL) || (p->contours == NULL))
2915     return FALSE;
2916   cur = p->contours;
2917   if (poly_InsideContour (cur, v0))
2918     {
2919       for (cur = cur->next; cur != NULL; cur = cur->next)
2920         if (poly_InsideContour (cur, v0))
2921           return FALSE;
2922       return TRUE;
2923     }
2924   return FALSE;
2925 }
2926
2927 BOOLp
2928 poly_M_CheckInside (POLYAREA * p, Vector v0)
2929 {
2930   POLYAREA *cur;
2931
2932   if ((p == NULL) || (v0 == NULL))
2933     return FALSE;
2934   cur = p;
2935   do
2936     {
2937       if (poly_CheckInside (cur, v0))
2938         return TRUE;
2939     }
2940   while ((cur = cur->f) != p);
2941   return FALSE;
2942 }
2943
2944 static double
2945 dot (Vector A, Vector B)
2946 {
2947   return (double)A[0] * (double)B[0] + (double)A[1] * (double)B[1];
2948 }
2949
2950 /* Compute whether point is inside a triangle formed by 3 other pints */
2951 /* Algorithm from http://www.blackpawn.com/texts/pointinpoly/default.html */
2952 static int
2953 point_in_triangle (Vector A, Vector B, Vector C, Vector P)
2954 {
2955   Vector v0, v1, v2;
2956   double dot00, dot01, dot02, dot11, dot12;
2957   double invDenom;
2958   double u, v;
2959
2960   /* Compute vectors */
2961   v0[0] = C[0] - A[0];  v0[1] = C[1] - A[1];
2962   v1[0] = B[0] - A[0];  v1[1] = B[1] - A[1];
2963   v2[0] = P[0] - A[0];  v2[1] = P[1] - A[1];
2964
2965   /* Compute dot products */
2966   dot00 = dot (v0, v0);
2967   dot01 = dot (v0, v1);
2968   dot02 = dot (v0, v2);
2969   dot11 = dot (v1, v1);
2970   dot12 = dot (v1, v2);
2971
2972   /* Compute barycentric coordinates */
2973   invDenom = 1. / (dot00 * dot11 - dot01 * dot01);
2974   u = (dot11 * dot02 - dot01 * dot12) * invDenom;
2975   v = (dot00 * dot12 - dot01 * dot02) * invDenom;
2976
2977   /* Check if point is in triangle */
2978   return (u > 0.0) && (v > 0.0) && (u + v < 1.0);
2979 }
2980
2981
2982 /* Returns the dot product of Vector A->B, and a vector
2983  * orthogonal to Vector C->D. The result is not normalisd, so will be
2984  * weighted by the magnitude of the C->D vector.
2985  */
2986 static double
2987 dot_orthogonal_to_direction (Vector A, Vector B, Vector C, Vector D)
2988 {
2989   Vector l1, l2, l3;
2990   l1[0] = B[0] - A[0];  l1[1] = B[1] - A[1];
2991   l2[0] = D[0] - C[0];  l2[1] = D[1] - C[1];
2992
2993   l3[0] = -l2[1];
2994   l3[1] = l2[0];
2995
2996   return dot (l1, l3);
2997 }
2998
2999 /* Algorithm from http://www.exaflop.org/docs/cgafaq/cga2.html
3000  *
3001  * "Given a simple polygon, find some point inside it. Here is a method based
3002  * on the proof that there exists an internal diagonal, in [O'Rourke, 13-14].
3003  * The idea is that the midpoint of a diagonal is interior to the polygon.
3004  *
3005  * 1.  Identify a convex vertex v; let its adjacent vertices be a and b.
3006  * 2.  For each other vertex q do:
3007  * 2a. If q is inside avb, compute distance to v (orthogonal to ab).
3008  * 2b. Save point q if distance is a new min.
3009  * 3.  If no point is inside, return midpoint of ab, or centroid of avb.
3010  * 4.  Else if some point inside, qv is internal: return its midpoint."
3011  *
3012  * [O'Rourke]: Computational Geometry in C (2nd Ed.)
3013  *             Joseph O'Rourke, Cambridge University Press 1998,
3014  *             ISBN 0-521-64010-5 Pbk, ISBN 0-521-64976-5 Hbk
3015  */
3016 static void
3017 poly_ComputeInteriorPoint (PLINE *poly, Vector v)
3018 {
3019   VNODE *pt1, *pt2, *pt3, *q;
3020   VNODE *min_q = NULL;
3021   double dist;
3022   double min_dist = 0.0;
3023   double dir = (poly->Flags.orient == PLF_DIR) ? 1. : -1;
3024
3025   /* Find a convex node on the polygon */
3026   pt1 = &poly->head;
3027   do
3028     {
3029       double dot_product;
3030
3031       pt2 = pt1->next;
3032       pt3 = pt2->next;
3033
3034       dot_product = dot_orthogonal_to_direction (pt1->point, pt2->point,
3035                                                  pt3->point, pt2->point);
3036
3037       if (dot_product * dir > 0.)
3038         break;
3039     }
3040   while ((pt1 = pt1->next) != &poly->head);
3041
3042   /* Loop over remaining vertices */
3043   q = pt3;
3044   do
3045     {
3046       /* Is current vertex "q" outside pt1 pt2 pt3 triangle? */
3047       if (!point_in_triangle (pt1->point, pt2->point, pt3->point, q->point))
3048         continue;
3049
3050       /* NO: compute distance to pt2 (v) orthogonal to pt1 - pt3) */
3051       /*     Record minimum */
3052       dist = dot_orthogonal_to_direction (q->point, pt2->point, pt1->point, pt3->point);
3053       if (min_q == NULL || dist < min_dist) {
3054         min_dist = dist;
3055         min_q = q;
3056       }
3057     }
3058   while ((q = q->next) != pt2);
3059
3060   /* Were any "q" found inside pt1 pt2 pt3? */
3061   if (min_q == NULL) {
3062     /* NOT FOUND: Return midpoint of pt1 pt3 */
3063     v[0] = (pt1->point[0] + pt3->point[0]) / 2;
3064     v[1] = (pt1->point[1] + pt3->point[1]) / 2;
3065   } else {
3066     /* FOUND: Return mid point of min_q, pt2 */
3067     v[0] = (pt2->point[0] + min_q->point[0]) / 2;
3068     v[1] = (pt2->point[1] + min_q->point[1]) / 2;
3069   }
3070 }
3071
3072
3073 /* NB: This function assumes the caller _knows_ the contours do not
3074  *     intersect. If the contours intersect, the result is undefined.
3075  *     It will return the correct result if the two contours share
3076  *     common points beteween their contours. (Identical contours
3077  *     are treated as being inside each other).
3078  */
3079 int
3080 poly_ContourInContour (PLINE * poly, PLINE * inner)
3081 {
3082   Vector point;
3083   assert (poly != NULL);
3084   assert (inner != NULL);
3085   if (cntrbox_inside (inner, poly))
3086     {
3087       /* We need to prove the "inner" contour is not outside
3088        * "poly" contour. If it is outside, we can return.
3089        */
3090       if (!poly_InsideContour (poly, inner->head.point))
3091         return 0;
3092
3093       poly_ComputeInteriorPoint (inner, point);
3094       return poly_InsideContour (poly, point);
3095     }
3096   return 0;
3097 }
3098
3099 void
3100 poly_Init (POLYAREA * p)
3101 {
3102   p->f = p->b = p;
3103   p->contours = NULL;
3104   p->contour_tree = r_create_tree (NULL, 0, 0);
3105 }
3106
3107 POLYAREA *
3108 poly_Create (void)
3109 {
3110   POLYAREA *res;
3111
3112   if ((res = (POLYAREA *)malloc (sizeof (POLYAREA))) != NULL)
3113     poly_Init (res);
3114   return res;
3115 }
3116
3117 void
3118 poly_FreeContours (PLINE ** pline)
3119 {
3120   PLINE *pl;
3121
3122   while ((pl = *pline) != NULL)
3123     {
3124       *pline = pl->next;
3125       poly_DelContour (&pl);
3126     }
3127 }
3128
3129 void
3130 poly_Free (POLYAREA ** p)
3131 {
3132   POLYAREA *cur;
3133
3134   if (*p == NULL)
3135     return;
3136   for (cur = (*p)->f; cur != *p; cur = (*p)->f)
3137     {
3138       poly_FreeContours (&cur->contours);
3139       r_destroy_tree (&cur->contour_tree);
3140       cur->f->b = cur->b;
3141       cur->b->f = cur->f;
3142       free (cur);
3143     }
3144   poly_FreeContours (&cur->contours);
3145   r_destroy_tree (&cur->contour_tree);
3146   free (*p), *p = NULL;
3147 }
3148
3149 static BOOLp
3150 inside_sector (VNODE * pn, Vector p2)
3151 {
3152   Vector cdir, ndir, pdir;
3153   int p_c, n_c, p_n;
3154
3155   assert (pn != NULL);
3156   vect_sub (cdir, p2, pn->point);
3157   vect_sub (pdir, pn->point, pn->prev->point);
3158   vect_sub (ndir, pn->next->point, pn->point);
3159
3160   p_c = vect_det2 (pdir, cdir) >= 0;
3161   n_c = vect_det2 (ndir, cdir) >= 0;
3162   p_n = vect_det2 (pdir, ndir) >= 0;
3163
3164   if ((p_n && p_c && n_c) || ((!p_n) && (p_c || n_c)))
3165     return TRUE;
3166   else
3167     return FALSE;
3168 }                               /* inside_sector */
3169
3170 /* returns TRUE if bad contour */
3171 BOOLp
3172 poly_ChkContour (PLINE * a)
3173 {
3174   VNODE *a1, *a2, *hit1, *hit2;
3175   Vector i1, i2;
3176   int icnt;
3177
3178   assert (a != NULL);
3179   a1 = &a->head;
3180   do
3181     {
3182       a2 = a1;
3183       do
3184         {
3185           if (!node_neighbours (a1, a2) &&
3186               (icnt = vect_inters2 (a1->point, a1->next->point,
3187                                     a2->point, a2->next->point, i1, i2)) > 0)
3188             {
3189               if (icnt > 1)
3190                 return TRUE;
3191
3192               if (vect_dist2 (i1, a1->point) < EPSILON)
3193                 hit1 = a1;
3194               else if (vect_dist2 (i1, a1->next->point) < EPSILON)
3195                 hit1 = a1->next;
3196               else
3197                 hit1 = NULL;
3198
3199               if (vect_dist2 (i1, a2->point) < EPSILON)
3200                 hit2 = a2;
3201               else if (vect_dist2 (i1, a2->next->point) < EPSILON)
3202                 hit2 = a2->next;
3203               else
3204                 hit2 = NULL;
3205
3206               if ((hit1 == NULL) && (hit2 == NULL))
3207                 {
3208                   /* If the intersection didn't land on an end-point of either
3209                    * line, we know the lines cross and we return TRUE.
3210                    */
3211                   return TRUE;
3212                 }
3213               else if (hit1 == NULL)
3214                 {
3215                 /* An end-point of the second line touched somewhere along the
3216                    length of the first line. Check where the second line leads. */
3217                   if (inside_sector (hit2, a1->point) !=
3218                       inside_sector (hit2, a1->next->point))
3219                     return TRUE;
3220                 }
3221               else if (hit2 == NULL)
3222                 {
3223                 /* An end-point of the first line touched somewhere along the
3224                    length of the second line. Check where the first line leads. */
3225                   if (inside_sector (hit1, a2->point) !=
3226                       inside_sector (hit1, a2->next->point))
3227                     return TRUE;
3228                 }
3229               else
3230                 {
3231                 /* Both lines intersect at an end-point. Check where they lead. */
3232                   if (inside_sector (hit1, hit2->prev->point) !=
3233                       inside_sector (hit1, hit2->next->point))
3234                     return TRUE;
3235                 }
3236             }
3237         }
3238       while ((a2 = a2->next) != &a->head);
3239     }
3240   while ((a1 = a1->next) != &a->head);
3241   return FALSE;
3242 }
3243
3244
3245 BOOLp
3246 poly_Valid (POLYAREA * p)
3247 {
3248   PLINE *c;
3249
3250   if ((p == NULL) || (p->contours == NULL))
3251     return FALSE;
3252
3253   if (p->contours->Flags.orient == PLF_INV || poly_ChkContour (p->contours))
3254     {
3255 #ifndef NDEBUG
3256       VNODE *v, *n;
3257       DEBUGP ("Invalid Outer PLINE\n");
3258       if (p->contours->Flags.orient == PLF_INV)
3259         DEBUGP ("failed orient\n");
3260       if (poly_ChkContour (p->contours))
3261         DEBUGP ("failed self-intersection\n");
3262       v = &p->contours->head;
3263       do
3264         {
3265           n = v->next;
3266           pcb_fprintf (stderr, "Line [%#mS %#mS %#mS %#mS 100 100 \"\"]\n",
3267                    v->point[0], v->point[1], n->point[0], n->point[1]);
3268         }
3269       while ((v = v->next) != &p->contours->head);
3270 #endif
3271       return FALSE;
3272     }
3273   for (c = p->contours->next; c != NULL; c = c->next)
3274     {
3275       if (c->Flags.orient == PLF_DIR ||
3276           poly_ChkContour (c) || !poly_ContourInContour (p->contours, c))
3277         {
3278 #ifndef NDEBUG
3279           VNODE *v, *n;
3280           DEBUGP ("Invalid Inner PLINE orient = %d\n", c->Flags.orient);
3281           if (c->Flags.orient == PLF_DIR)
3282             DEBUGP ("failed orient\n");
3283           if (poly_ChkContour (c))
3284             DEBUGP ("failed self-intersection\n");
3285           if (!poly_ContourInContour (p->contours, c))
3286             DEBUGP ("failed containment\n");
3287           v = &c->head;
3288           do
3289             {
3290               n = v->next;
3291               pcb_fprintf (stderr, "Line [%#mS %#mS %#mS %#mS 100 100 \"\"]\n",
3292                        v->point[0], v->point[1], n->point[0], n->point[1]);
3293             }
3294           while ((v = v->next) != &c->head);
3295 #endif
3296           return FALSE;
3297         }
3298     }
3299   return TRUE;
3300 }
3301
3302
3303 Vector vect_zero = { (long) 0, (long) 0 };
3304
3305 /*********************************************************************/
3306 /*             L o n g   V e c t o r   S t u f f                     */
3307 /*********************************************************************/
3308
3309 void
3310 vect_init (Vector v, double x, double y)
3311 {
3312   v[0] = (long) x;
3313   v[1] = (long) y;
3314 }                               /* vect_init */
3315
3316 #define Vzero(a)   ((a)[0] == 0. && (a)[1] == 0.)
3317
3318 #define Vsub(a,b,c) {(a)[0]=(b)[0]-(c)[0];(a)[1]=(b)[1]-(c)[1];}
3319
3320 int
3321 vect_equal (Vector v1, Vector v2)
3322 {
3323   return (v1[0] == v2[0] && v1[1] == v2[1]);
3324 }                               /* vect_equal */
3325
3326
3327 void
3328 vect_sub (Vector res, Vector v1, Vector v2)
3329 {
3330 Vsub (res, v1, v2)}             /* vect_sub */
3331
3332 void
3333 vect_min (Vector v1, Vector v2, Vector v3)
3334 {
3335   v1[0] = (v2[0] < v3[0]) ? v2[0] : v3[0];
3336   v1[1] = (v2[1] < v3[1]) ? v2[1] : v3[1];
3337 }                               /* vect_min */
3338
3339 void
3340 vect_max (Vector v1, Vector v2, Vector v3)
3341 {
3342   v1[0] = (v2[0] > v3[0]) ? v2[0] : v3[0];
3343   v1[1] = (v2[1] > v3[1]) ? v2[1] : v3[1];
3344 }                               /* vect_max */
3345
3346 double
3347 vect_len2 (Vector v)
3348 {
3349   return ((double) v[0] * v[0] + (double) v[1] * v[1]); /* why sqrt? only used for compares */
3350 }
3351
3352 double
3353 vect_dist2 (Vector v1, Vector v2)
3354 {
3355   double dx = v1[0] - v2[0];
3356   double dy = v1[1] - v2[1];
3357
3358   return (dx * dx + dy * dy);   /* why sqrt */
3359 }
3360
3361 /* value has sign of angle between vectors */
3362 double
3363 vect_det2 (Vector v1, Vector v2)
3364 {
3365   return (((double) v1[0] * v2[1]) - ((double) v2[0] * v1[1]));
3366 }
3367
3368 static double
3369 vect_m_dist (Vector v1, Vector v2)
3370 {
3371   double dx = v1[0] - v2[0];
3372   double dy = v1[1] - v2[1];
3373   double dd = (dx * dx + dy * dy);      /* sqrt */
3374
3375   if (dx > 0)
3376     return +dd;
3377   if (dx < 0)
3378     return -dd;
3379   if (dy > 0)
3380     return +dd;
3381   return -dd;
3382 }                               /* vect_m_dist */
3383
3384 /*
3385 vect_inters2
3386  (C) 1993 Klamer Schutte
3387  (C) 1997 Michael Leonov, Alexey Nikitin
3388 */
3389
3390 int
3391 vect_inters2 (Vector p1, Vector p2, Vector q1, Vector q2,
3392               Vector S1, Vector S2)
3393 {
3394   double s, t, deel;
3395   double rpx, rpy, rqx, rqy;
3396
3397   if (max (p1[0], p2[0]) < min (q1[0], q2[0]) ||
3398       max (q1[0], q2[0]) < min (p1[0], p2[0]) ||
3399       max (p1[1], p2[1]) < min (q1[1], q2[1]) ||
3400       max (q1[1], q2[1]) < min (p1[1], p2[1]))
3401     return 0;
3402
3403   rpx = p2[0] - p1[0];
3404   rpy = p2[1] - p1[1];
3405   rqx = q2[0] - q1[0];
3406   rqy = q2[1] - q1[1];
3407
3408   deel = rpy * rqx - rpx * rqy; /* -vect_det(rp,rq); */
3409
3410   /* coordinates are 30-bit integers so deel will be exactly zero
3411    * if the lines are parallel
3412    */
3413
3414   if (deel == 0)                /* parallel */
3415     {
3416       double dc1, dc2, d1, d2, h;       /* Check to see whether p1-p2 and q1-q2 are on the same line */
3417       Vector hp1, hq1, hp2, hq2, q1p1, q1q2;
3418
3419       Vsub2 (q1p1, q1, p1);
3420       Vsub2 (q1q2, q1, q2);
3421
3422
3423       /* If this product is not zero then p1-p2 and q1-q2 are not on same line! */
3424       if (vect_det2 (q1p1, q1q2) != 0)
3425         return 0;
3426       dc1 = 0;                  /* m_len(p1 - p1) */
3427
3428       dc2 = vect_m_dist (p1, p2);
3429       d1 = vect_m_dist (p1, q1);
3430       d2 = vect_m_dist (p1, q2);
3431
3432 /* Sorting the independent points from small to large */
3433       Vcpy2 (hp1, p1);
3434       Vcpy2 (hp2, p2);
3435       Vcpy2 (hq1, q1);
3436       Vcpy2 (hq2, q2);
3437       if (dc1 > dc2)
3438         {                       /* hv and h are used as help-variable. */
3439           Vswp2 (hp1, hp2);
3440           h = dc1, dc1 = dc2, dc2 = h;
3441         }
3442       if (d1 > d2)
3443         {
3444           Vswp2 (hq1, hq2);
3445           h = d1, d1 = d2, d2 = h;
3446         }
3447
3448 /* Now the line-pieces are compared */
3449
3450       if (dc1 < d1)
3451         {
3452           if (dc2 < d1)
3453             return 0;
3454           if (dc2 < d2)
3455             {
3456               Vcpy2 (S1, hp2);
3457               Vcpy2 (S2, hq1);
3458             }
3459           else
3460             {
3461               Vcpy2 (S1, hq1);
3462               Vcpy2 (S2, hq2);
3463             };
3464         }
3465       else
3466         {
3467           if (dc1 > d2)
3468             return 0;
3469           if (dc2 < d2)
3470             {
3471               Vcpy2 (S1, hp1);
3472               Vcpy2 (S2, hp2);
3473             }
3474           else
3475             {
3476               Vcpy2 (S1, hp1);
3477               Vcpy2 (S2, hq2);
3478             };
3479         }
3480       return (Vequ2 (S1, S2) ? 1 : 2);
3481     }
3482   else
3483     {                           /* not parallel */
3484       /*
3485        * We have the lines:
3486        * l1: p1 + s(p2 - p1)
3487        * l2: q1 + t(q2 - q1)
3488        * And we want to know the intersection point.
3489        * Calculate t:
3490        * p1 + s(p2-p1) = q1 + t(q2-q1)
3491        * which is similar to the two equations:
3492        * p1x + s * rpx = q1x + t * rqx
3493        * p1y + s * rpy = q1y + t * rqy
3494        * Multiplying these by rpy resp. rpx gives:
3495        * rpy * p1x + s * rpx * rpy = rpy * q1x + t * rpy * rqx
3496        * rpx * p1y + s * rpx * rpy = rpx * q1y + t * rpx * rqy
3497        * Subtracting these gives:
3498        * rpy * p1x - rpx * p1y = rpy * q1x - rpx * q1y + t * ( rpy * rqx - rpx * rqy )
3499        * So t can be isolated:
3500        * t = (rpy * ( p1x - q1x ) + rpx * ( - p1y + q1y )) / ( rpy * rqx - rpx * rqy )
3501        * and s can be found similarly
3502        * s = (rqy * (q1x - p1x) + rqx * (p1y - q1y))/( rqy * rpx - rqx * rpy)
3503        */
3504
3505       if (Vequ2 (q1, p1) || Vequ2 (q1, p2))
3506         {
3507           S1[0] = q1[0];
3508           S1[1] = q1[1];
3509         }
3510       else if (Vequ2 (q2, p1) || Vequ2 (q2, p2))
3511         {
3512           S1[0] = q2[0];
3513           S1[1] = q2[1];
3514         }
3515       else
3516         {
3517           s = (rqy * (p1[0] - q1[0]) + rqx * (q1[1] - p1[1])) / deel;
3518           if (s < 0 || s > 1.)
3519             return 0;
3520           t = (rpy * (p1[0] - q1[0]) + rpx * (q1[1] - p1[1])) / deel;
3521           if (t < 0 || t > 1.)
3522             return 0;
3523
3524           S1[0] = q1[0] + ROUND (t * rqx);
3525           S1[1] = q1[1] + ROUND (t * rqy);
3526         }
3527       return 1;
3528     }
3529 }                               /* vect_inters2 */
3530
3531 /* how about expanding polygons so that edges can be arcs rather than
3532  * lines. Consider using the third coordinate to store the radius of the
3533  * arc. The arc would pass through the vertex points. Positive radius
3534  * would indicate the arc bows left (center on right of P1-P2 path)
3535  * negative radius would put the center on the other side. 0 radius
3536  * would mean the edge is a line instead of arc
3537  * The intersections of the two circles centered at the vertex points
3538  * would determine the two possible arc centers. If P2.x > P1.x then
3539  * the center with smaller Y is selected for positive r. If P2.y > P1.y
3540  * then the center with greate X is selected for positive r.
3541  *
3542  * the vec_inters2() routine would then need to handle line-line
3543  * line-arc and arc-arc intersections.
3544  *
3545  * perhaps reverse tracing the arc would require look-ahead to check
3546  * for arcs
3547  */