src/polygon1.c

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