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