gpc.d

   1 /**
   2 Generic Polygon Clipper
   3
   4 A new algorithm for calculating the difference, intersection,
   5 exclusive-or or union of arbitrary polygon sets.
   6
   7 Written by Alan Murta (email: gpc@cs.man.ac.uk)
   8
   9 [http://www.cs.man.ac.uk/~toby/alan/software/|Original GPC site].
  10
  11 Current version is 2.33 (21st May 2014).
  12
  13 Ported to D by Ketmar // Invisible Vector.
  14
  15 GPC is using standard libc `malloc()`, `realloc()` and `free()` functions
  16 to manage memory.
  17
  18 Copyright (C) Advanced Interfaces Group, University of Manchester.
  19
  20 This software is free for non-commercial use. It may be copied,
  21 modified, and redistributed provided that this copyright notice
  22 is preserved on all copies. The intellectual property rights of
  23 the algorithms used reside with the University of Manchester
  24 Advanced Interfaces Group.
  25
  26 You may not use this software, in whole or in part, in support
  27 of any commercial product without the express consent of the
  28 author.
  29
  30 There is no warranty or other guarantee of fitness of this
  31 software for any purpose. It is provided solely "as is".
  32 */
  33 module iv.gpc;
  34 private nothrow @trusted @nogc:
  35
  36 /*
  37 ===========================================================================
  38                                Constants
  39 ===========================================================================
  40 */
  41
  42 private enum DBL_EPSILON = double.epsilon; //cast(double)2.2204460492503131e-16;
  43
  44 version(gpc_use_custom_epsilon) {
  45   public double GPC_EPSILON = DBL_EPSILON; /// Increase GPC_EPSILON to encourage merging of near coincident edges
  46 } else {
  47   public enum GPC_EPSILON = DBL_EPSILON; /// Increase GPC_EPSILON to encourage merging of near coincident edges
  48 }
  49
  50 public bool GPC_INVERT_TRISTRIPS = false;
  51
  52 version(gpc_use_custom_epsilon) {
  53   /// Compares two floating point numbers using [GPC_EPSILON] as epsilon value.
  54   public bool gpc_eq (in double a, in double b) nothrow @safe @nogc { pragma(inline, true); import std.math; return (abs(a-b) <= GPC_EPSILON); }
  55 } else {
  56   /// Compares two floating point numbers using [GPC_EPSILON] as epsilon value.
  57   public bool gpc_eq (in double a, in double b) pure nothrow @safe @nogc { pragma(inline, true); import std.math; return (abs(a-b) <= GPC_EPSILON); }
  58 }
  59
  60
  61 /*
  62 ===========================================================================
  63                            Public Data Types
  64 ===========================================================================
  65 */
  66
  67 /// Set operation type
  68 public enum GPC {
  69   Diff, /// Difference
  70   Int, /// Intersection
  71   Xor, /// Exclusive or
  72   Union, /// Union
  73 }
  74
  75 /// Polygon vertex structure
  76 public struct gpc_vertex {
  77   double x; /// vertex x component
  78   double y; /// vertex y component
  79 }
  80
  81 /// Vertex list structure
  82 public struct gpc_vertex_list {
  83   int num_vertices; /// Number of vertices in list
  84   gpc_vertex* vertex; /// Vertex array pointer
  85 }
  86
  87 /// Polygon set structure
  88 public struct gpc_polygon {
  89   int num_contours; /// Number of contours in polygon
  90   int* hole; /// Hole / external contour flags
  91   gpc_vertex_list* contour; /// Contour array pointer
  92 }
  93
  94 /// Tristrip set structure
  95 public struct gpc_tristrip {
  96   int num_strips; /// Number of tristrips
  97   gpc_vertex_list* strip; /// Tristrip array pointer
  98 }
  99
 100
 101 private:
 102 /*
 103 ===========================================================================
 104                                 Constants
 105 ===========================================================================
 106 */
 107
 108 enum LEFT = 0;
 109 enum RIGHT = 1;
 110
 111 enum ABOVE = 0;
 112 enum BELOW = 1;
 113
 114 enum CLIP = 0;
 115 enum SUBJ = 1;
 116
 117
 118 /*
 119 ===========================================================================
 120                                  Macros
 121 ===========================================================================
 122 */
 123
 124 version(gpc_use_custom_epsilon) {
 125   bool EQ(T) (in T a, in T b) nothrow @safe @nogc { pragma(inline, true); import std.math; return (abs(a-b) <= GPC_EPSILON); }
 126 } else {
 127   bool EQ(T) (in T a, in T b) pure nothrow @safe @nogc { pragma(inline, true); import std.math; return (abs(a-b) <= GPC_EPSILON); }
 128 }
 129
 130 T PREV_INDEX(T) (in T i, in T n) pure nothrow @safe @nogc { pragma(inline, true); return (i-1+n)%n; }
 131 T NEXT_INDEX(T) (in T i, in T n) pure nothrow @safe @nogc { pragma(inline, true); return (i+1)%n; }
 132
 133 bool OPTIMAL(TV, T) (const(TV)* v, in T i, in T n) pure nothrow @trusted @nogc { pragma(inline, true); return ((v[PREV_INDEX(i, n)].y != v[i].y) || (v[NEXT_INDEX(i, n)].y != v[i].y)); }
 134
 135 bool FWD_MIN(TV, T) (const(TV)* v, in T i, in T n) pure nothrow @trusted @nogc { pragma(inline, true); return ((v[PREV_INDEX(i, n)].vertex.y >= v[i].vertex.y) && (v[NEXT_INDEX(i, n)].vertex.y > v[i].vertex.y)); }
 136
 137 bool NOT_FMAX(TV, T) (const(TV)* v, in T i, in T n) pure nothrow @trusted @nogc { pragma(inline, true); return (v[NEXT_INDEX(i, n)].vertex.y > v[i].vertex.y); }
 138
 139 bool REV_MIN(TV, T) (const(TV)* v, in T i, in T n) pure nothrow @trusted @nogc { pragma(inline, true); return ((v[PREV_INDEX(i, n)].vertex.y > v[i].vertex.y) && (v[NEXT_INDEX(i, n)].vertex.y >= v[i].vertex.y)); }
 140
 141 bool NOT_RMAX(TV, T) (const(TV)* v, in T i, in T n) pure nothrow @trusted @nogc { pragma(inline, true); return (v[PREV_INDEX(i, n)].vertex.y > v[i].vertex.y); }
 142
 143 void VERTEX(T) (ref edge_node* e, in int p, in int s, in T x, in T y) nothrow @trusted @nogc { pragma(inline, true); add_vertex(&((e).outp[(p)].v[(s)]), x, y); (e).outp[(p)].active++; }
 144
 145 void P_EDGE(T) (ref edge_node* d, ref edge_node* e, in int p, ref T i, in T j) nothrow @trusted @nogc {
 146   (d)= (e);
 147   do { (d)= (d).prev; } while (!(d).outp[(p)]);
 148   (i)= (d).bot.x + (d).dx * ((j)-(d).bot.y);
 149 }
 150
 151 void N_EDGE(T) (ref edge_node* d, ref edge_node* e, in int p, ref T i, in T j) nothrow @trusted @nogc {
 152   (d) = (e);
 153   do { (d)= (d).next; } while (!(d).outp[(p)]);
 154   (i) = (d).bot.x + (d).dx * ((j)-(d).bot.y);
 155 }
 156
 157
 158 T* MALLOC(T) (uint count=1) nothrow @trusted @nogc {
 159   import core.stdc.stdlib : malloc;
 160   static assert(T.sizeof > 0);
 161   if (count == 0) return null;
 162   if (count >= int.max/4/T.sizeof) assert(0, "gpc malloc failure");
 163   T* p = cast(T*)malloc(T.sizeof*count);
 164   if (p is null) assert(0, "gpc malloc failure");
 165   return p;
 166 }
 167
 168 void REALLOC(T) (ref T* p, uint count) nothrow @trusted @nogc {
 169   import core.stdc.stdlib : free, realloc;
 170   if (p is null) {
 171     if (count != 0) p = MALLOC!T(count);
 172   } else if (count == 0) {
 173     if (p !is null) {
 174       free(p);
 175       p = null;
 176     }
 177   } else {
 178     static assert(T.sizeof > 0);
 179     assert(count > 0);
 180     if (count >= int.max/4/T.sizeof) assert(0, "gpc malloc failure");
 181     T* np = cast(T*)realloc(p, T.sizeof*count);
 182     if (np is null) assert(0, "gpc malloc failure");
 183     p = np;
 184   }
 185 }
 186
 187 void FREE(T) (ref T* ptr) nothrow @trusted @nogc {
 188   import core.stdc.stdlib : free;
 189   if (ptr !is null) { free(ptr); ptr = null; }
 190 }
 191
 192
 193 /*
 194 ===========================================================================
 195                             Private Data Types
 196 ===========================================================================
 197 */
 198
 199 // Edge intersection classes
 200 alias vertex_type = int;
 201 enum : uint {
 202   NUL, // Empty non-intersection
 203   EMX, // External maximum
 204   ELI, // External left intermediate
 205   TED, // Top edge
 206   ERI, // External right intermediate
 207   RED, // Right edge
 208   IMM, // Internal maximum and minimum
 209   IMN, // Internal minimum
 210   EMN, // External minimum
 211   EMM, // External maximum and minimum
 212   LED, // Left edge
 213   ILI, // Internal left intermediate
 214   BED, // Bottom edge
 215   IRI, // Internal right intermediate
 216   IMX, // Internal maximum
 217   FUL, // Full non-intersection
 218 }
 219
 220 // Horizontal edge states
 221 alias h_state = int;
 222 enum : int {
 223   NH, // No horizontal edge
 224   BH, // Bottom horizontal edge
 225   TH, // Top horizontal edge
 226 }
 227
 228 // Edge bundle state
 229 alias bundle_state = int;
 230 enum : int {
 231   UNBUNDLED, // Isolated edge not within a bundle
 232   BUNDLE_HEAD, // Bundle head node
 233   BUNDLE_TAIL, // Passive bundle tail node
 234 }
 235
 236 // Internal vertex list datatype
 237 //alias v_shape = vertex_node;
 238 struct vertex_node {
 239   double x; // X coordinate component
 240   double y; // Y coordinate component
 241   vertex_node* next; // Pointer to next vertex in list
 242 }
 243
 244 //p_shape
 245 // Internal contour / tristrip type
 246 struct polygon_node {
 247   int active; // Active flag / vertex count
 248   int hole; // Hole / external contour flag
 249   vertex_node*[2] v; // Left and right vertex list ptrs
 250   polygon_node* next; // Pointer to next polygon contour
 251   polygon_node* proxy; // Pointer to actual structure used
 252 }
 253
 254 // edge_shape
 255 struct edge_node {
 256   gpc_vertex vertex; // Piggy-backed contour vertex data
 257   gpc_vertex bot; // Edge lower (x, y) coordinate
 258   gpc_vertex top; // Edge upper (x, y) coordinate
 259   double xb; // Scanbeam bottom x coordinate
 260   double xt; // Scanbeam top x coordinate
 261   double dx; // Change in x for a unit y increase
 262   int type; // Clip / subject edge flag
 263   int[2][2] bundle; // Bundle edge flags
 264   int[2] bside; // Bundle left / right indicators
 265   bundle_state[2] bstate; // Edge bundle state
 266   polygon_node*[2] outp; // Output polygon / tristrip pointer
 267   edge_node* prev; // Previous edge in the AET
 268   edge_node* next; // Next edge in the AET
 269   edge_node* pred; // Edge connected at the lower end
 270   edge_node* succ; // Edge connected at the upper end
 271   edge_node* next_bound; // Pointer to next bound in LMT
 272 }
 273
 274 // lmt_shape
 275 // Local minima table
 276 struct lmt_node {
 277   double y; // Y coordinate at local minimum
 278   edge_node* first_bound; // Pointer to bound list
 279   lmt_node* next; // Pointer to next local minimum
 280 }
 281
 282 // sbt_t_shape
 283 // Scanbeam tree
 284 struct sb_tree {
 285   double y; // Scanbeam node y value
 286   sb_tree* less; // Pointer to nodes with lower y
 287   sb_tree* more; // Pointer to nodes with higher y
 288 }
 289
 290 // it_shape
 291 // Intersection table
 292 struct it_node {
 293   edge_node*[2] ie; // Intersecting edge (bundle) pair
 294   gpc_vertex point; // Point of intersection
 295   it_node* next; // The next intersection table node
 296 }
 297
 298 // st_shape
 299 // Sorted edge table
 300 struct st_node {
 301   edge_node* edge; // Pointer to AET edge
 302   double xb; // Scanbeam bottom x coordinate
 303   double xt; // Scanbeam top x coordinate
 304   double dx; // Change in x for a unit y increase
 305   st_node* prev; // Previous edge in sorted list
 306 }
 307
 308 // bbox_shape
 309 // Contour axis-aligned bounding box
 310 struct bbox {
 311   double xmin; // Minimum x coordinate
 312   double ymin; // Minimum y coordinate
 313   double xmax; // Maximum x coordinate
 314   double ymax; // Maximum y coordinate
 315 }
 316
 317
 318 /*
 319 ===========================================================================
 320                                Global Data
 321 ===========================================================================
 322 */
 323
 324 // Horizontal edge state transitions within scanbeam boundary
 325 static immutable h_state[6][3] next_h_state= [
 326   /*        ABOVE     BELOW     CROSS */
 327   /*        L   R     L   R     L   R */
 328   /* NH */ [BH, TH,   TH, BH,   NH, NH],
 329   /* BH */ [NH, NH,   NH, NH,   TH, TH],
 330   /* TH */ [NH, NH,   NH, NH,   BH, BH],
 331 ];
 332
 333
 334 /*
 335 ===========================================================================
 336                              Private Functions
 337 ===========================================================================
 338 */
 339
 340 private void reset_it (it_node** it) {
 341   while (*it) {
 342     it_node* itn = (*it).next;
 343     FREE(*it);
 344     *it = itn;
 345   }
 346 }
 347
 348
 349 private void reset_lmt (lmt_node** lmt) {
 350   while (*lmt) {
 351     lmt_node* lmtn = (*lmt).next;
 352     FREE(*lmt);
 353     *lmt= lmtn;
 354   }
 355 }
 356
 357
 358 private void insert_bound (edge_node** b, edge_node* e) {
 359   edge_node* existing_bound;
 360
 361   if (!*b) {
 362     // Link node e to the tail of the list
 363     *b = e;
 364   } else {
 365     // Do primary sort on the x field
 366     if (e[0].bot.x < (*b)[0].bot.x) {
 367       // Insert a new node mid-list
 368       existing_bound = *b;
 369       *b = e;
 370       (*b).next_bound = existing_bound;
 371     } else {
 372       if (e[0].bot.x == (*b)[0].bot.x) {
 373         // Do secondary sort on the dx field
 374         if (e[0].dx < (*b)[0].dx) {
 375           // Insert a new node mid-list
 376           existing_bound = *b;
 377           *b = e;
 378           (*b).next_bound = existing_bound;
 379         } else {
 380           // Head further down the list
 381           insert_bound(&((*b).next_bound), e);
 382         }
 383       } else {
 384         // Head further down the list
 385         insert_bound(&((*b).next_bound), e);
 386       }
 387     }
 388   }
 389 }
 390
 391
 392 private edge_node** bound_list (lmt_node** lmt, double y) {
 393   lmt_node* existing_node;
 394
 395   if (!*lmt) {
 396     // Add node onto the tail end of the LMT
 397     *lmt = MALLOC!lmt_node();
 398     (*lmt).y = y;
 399     (*lmt).first_bound = null;
 400     (*lmt).next = null;
 401     return &((*lmt).first_bound);
 402   } else if (y < (*lmt).y) {
 403     // Insert a new LMT node before the current node
 404     existing_node = *lmt;
 405     *lmt = MALLOC!lmt_node();
 406     (*lmt).y = y;
 407     (*lmt).first_bound = null;
 408     (*lmt).next = existing_node;
 409     return &((*lmt).first_bound);
 410   } else if (y > (*lmt).y) {
 411     // Head further up the LMT
 412     return bound_list(&((*lmt).next), y);
 413   } else {
 414     // Use this existing LMT node
 415     return &((*lmt).first_bound);
 416   }
 417 }
 418
 419
 420 private void add_to_sbtree (int* entries, sb_tree** sbtree, double y) {
 421   if (!*sbtree) {
 422     // Add a new tree node here
 423     *sbtree = MALLOC!sb_tree();
 424     (*sbtree).y = y;
 425     (*sbtree).less = null;
 426     (*sbtree).more = null;
 427     ++(*entries);
 428   } else {
 429     if ((*sbtree).y > y) {
 430       // Head into the 'less' sub-tree
 431       add_to_sbtree(entries, &((*sbtree).less), y);
 432     } else {
 433       if ((*sbtree).y < y) {
 434         // Head into the 'more' sub-tree
 435         add_to_sbtree(entries, &((*sbtree).more), y);
 436       }
 437     }
 438   }
 439 }
 440
 441
 442 private void build_sbt (int* entries, double* sbt, sb_tree* sbtree) {
 443   if (sbtree.less) build_sbt(entries, sbt, sbtree.less);
 444   sbt[*entries]= sbtree.y;
 445   ++(*entries);
 446   if (sbtree.more) build_sbt(entries, sbt, sbtree.more);
 447 }
 448
 449
 450 private void free_sbtree (sb_tree** sbtree) {
 451   if (*sbtree) {
 452     free_sbtree(&((*sbtree).less));
 453     free_sbtree(&((*sbtree).more));
 454     FREE(*sbtree);
 455   }
 456 }
 457
 458
 459 private int count_optimal_vertices (gpc_vertex_list c) {
 460   int result = 0;
 461   // Ignore non-contributing contours
 462   if (c.num_vertices > 0) {
 463     foreach (immutable int i; 0..c.num_vertices) {
 464       // Ignore superfluous vertices embedded in horizontal edges
 465       if (OPTIMAL(c.vertex, i, c.num_vertices)) ++result;
 466     }
 467   }
 468   return result;
 469 }
 470
 471
 472 private edge_node* build_lmt (lmt_node** lmt, sb_tree** sbtree, int* sbt_entries, gpc_polygon* p, int type, GPC op) {
 473   int min, max, num_edges, v, num_vertices;
 474   int total_vertices = 0, e_index = 0;
 475   edge_node* e, edge_table;
 476
 477   foreach (immutable int c; 0..p.num_contours) total_vertices += count_optimal_vertices(p.contour[c]);
 478
 479   // Create the entire input polygon edge table in one go
 480   edge_table = MALLOC!edge_node(total_vertices);
 481
 482   foreach (immutable int c; 0..p.num_contours) {
 483     if (p.contour[c].num_vertices < 0) {
 484       // Ignore the non-contributing contour and repair the vertex count
 485       p.contour[c].num_vertices = -p.contour[c].num_vertices;
 486     } else {
 487       // Perform contour optimisation
 488       num_vertices = 0;
 489       foreach (immutable int i; 0..p.contour[c].num_vertices) {
 490         if (OPTIMAL(p.contour[c].vertex, i, p.contour[c].num_vertices)) {
 491           edge_table[num_vertices].vertex.x = p.contour[c].vertex[i].x;
 492           edge_table[num_vertices].vertex.y = p.contour[c].vertex[i].y;
 493           // Record vertex in the scanbeam table
 494           add_to_sbtree(sbt_entries, sbtree, edge_table[num_vertices].vertex.y);
 495           ++num_vertices;
 496         }
 497       }
 498
 499       // Do the contour forward pass
 500       for (min = 0; min < num_vertices; ++min) {
 501         // If a forward local minimum...
 502         if (FWD_MIN(edge_table, min, num_vertices)) {
 503           // Search for the next local maximum...
 504           num_edges = 1;
 505           max = NEXT_INDEX(min, num_vertices);
 506           while (NOT_FMAX(edge_table, max, num_vertices)) {
 507             ++num_edges;
 508             max = NEXT_INDEX(max, num_vertices);
 509           }
 510
 511           // Build the next edge list
 512           e = &edge_table[e_index];
 513           e_index += num_edges;
 514           v = min;
 515           e[0].bstate[BELOW] = UNBUNDLED;
 516           e[0].bundle[BELOW][CLIP] = false;
 517           e[0].bundle[BELOW][SUBJ] = false;
 518           foreach (immutable int i; 0..num_edges) {
 519             e[i].xb = edge_table[v].vertex.x;
 520             e[i].bot.x = edge_table[v].vertex.x;
 521             e[i].bot.y = edge_table[v].vertex.y;
 522
 523             v = NEXT_INDEX(v, num_vertices);
 524
 525             e[i].top.x = edge_table[v].vertex.x;
 526             e[i].top.y = edge_table[v].vertex.y;
 527             e[i].dx = (edge_table[v].vertex.x - e[i].bot.x) / (e[i].top.y - e[i].bot.y);
 528             e[i].type = type;
 529             e[i].outp[ABOVE] = null;
 530             e[i].outp[BELOW] = null;
 531             e[i].next = null;
 532             e[i].prev = null;
 533             e[i].succ = ((num_edges > 1) && (i < (num_edges - 1))) ? &(e[i + 1]) : null;
 534             e[i].pred = ((num_edges > 1) && (i > 0)) ? &(e[i - 1]) : null;
 535             e[i].next_bound = null;
 536             e[i].bside[CLIP] = (op == GPC.Diff ? RIGHT : LEFT);
 537             e[i].bside[SUBJ] = LEFT;
 538           }
 539           insert_bound(bound_list(lmt, edge_table[min].vertex.y), e);
 540         }
 541       }
 542
 543       // Do the contour reverse pass
 544       for (min = 0; min < num_vertices; ++min) {
 545         // If a reverse local minimum...
 546         if (REV_MIN(edge_table, min, num_vertices)) {
 547           // Search for the previous local maximum...
 548           num_edges = 1;
 549           max = PREV_INDEX(min, num_vertices);
 550           while (NOT_RMAX(edge_table, max, num_vertices)) {
 551             ++num_edges;
 552             max = PREV_INDEX(max, num_vertices);
 553           }
 554
 555           // Build the previous edge list
 556           e = &edge_table[e_index];
 557           e_index += num_edges;
 558           v = min;
 559           e[0].bstate[BELOW] = UNBUNDLED;
 560           e[0].bundle[BELOW][CLIP] = false;
 561           e[0].bundle[BELOW][SUBJ] = false;
 562           foreach (immutable i; 0..num_edges) {
 563             e[i].xb = edge_table[v].vertex.x;
 564             e[i].bot.x = edge_table[v].vertex.x;
 565             e[i].bot.y = edge_table[v].vertex.y;
 566
 567             v = PREV_INDEX(v, num_vertices);
 568
 569             e[i].top.x = edge_table[v].vertex.x;
 570             e[i].top.y = edge_table[v].vertex.y;
 571             e[i].dx = (edge_table[v].vertex.x - e[i].bot.x) / (e[i].top.y - e[i].bot.y);
 572             e[i].type = type;
 573             e[i].outp[ABOVE] = null;
 574             e[i].outp[BELOW] = null;
 575             e[i].next = null;
 576             e[i].prev = null;
 577             e[i].succ = ((num_edges > 1) && (i < (num_edges - 1))) ? &(e[i + 1]) : null;
 578             e[i].pred = ((num_edges > 1) && (i > 0)) ? &(e[i - 1]) : null;
 579             e[i].next_bound = null;
 580             e[i].bside[CLIP] = (op == GPC.Diff ? RIGHT : LEFT);
 581             e[i].bside[SUBJ] = LEFT;
 582           }
 583           insert_bound(bound_list(lmt, edge_table[min].vertex.y), e);
 584         }
 585       }
 586     }
 587   }
 588   return edge_table;
 589 }
 590
 591
 592 private void add_edge_to_aet (edge_node** aet, edge_node* edge, edge_node* prev) {
 593   if (!*aet) {
 594     // Append edge onto the tail end of the AET
 595     *aet = edge;
 596     edge.prev = prev;
 597     edge.next = null;
 598   } else {
 599     // Do primary sort on the xb field
 600     if (edge.xb < (*aet).xb) {
 601       // Insert edge here (before the AET edge)
 602       edge.prev = prev;
 603       edge.next = *aet;
 604       (*aet).prev = edge;
 605       *aet = edge;
 606     } else {
 607       if (edge.xb == (*aet).xb) {
 608         // Do secondary sort on the dx field
 609         if (edge.dx < (*aet).dx) {
 610           // Insert edge here (before the AET edge)
 611           edge.prev = prev;
 612           edge.next = *aet;
 613           (*aet).prev = edge;
 614           *aet = edge;
 615         } else {
 616           // Head further into the AET
 617           add_edge_to_aet(&((*aet).next), edge, *aet);
 618         }
 619       } else {
 620         // Head further into the AET
 621         add_edge_to_aet(&((*aet).next), edge, *aet);
 622       }
 623     }
 624   }
 625 }
 626
 627
 628 private void add_intersection (it_node** it, edge_node* edge0, edge_node* edge1, double x, double y) {
 629   it_node* existing_node;
 630
 631   if (!*it) {
 632     // Append a new node to the tail of the list
 633     *it = MALLOC!it_node();
 634     (*it).ie[0] = edge0;
 635     (*it).ie[1] = edge1;
 636     (*it).point.x = x;
 637     (*it).point.y = y;
 638     (*it).next = null;
 639   } else {
 640     if ((*it).point.y > y) {
 641       // Insert a new node mid-list
 642       existing_node = *it;
 643       *it = MALLOC!it_node();
 644       (*it).ie[0] = edge0;
 645       (*it).ie[1] = edge1;
 646       (*it).point.x = x;
 647       (*it).point.y = y;
 648       (*it).next = existing_node;
 649     } else {
 650       // Head further down the list
 651       add_intersection(&((*it).next), edge0, edge1, x, y);
 652     }
 653   }
 654 }
 655
 656
 657 private void add_st_edge (st_node** st, it_node** it, edge_node* edge, double dy) {
 658   import std.math : abs;
 659
 660   st_node* existing_node;
 661   double den, r, x, y;
 662
 663   if (!*st) {
 664     // Append edge onto the tail end of the ST
 665     *st = MALLOC!st_node();
 666     (*st).edge = edge;
 667     (*st).xb = edge.xb;
 668     (*st).xt = edge.xt;
 669     (*st).dx = edge.dx;
 670     (*st).prev = null;
 671   } else {
 672     den = ((*st).xt - (*st).xb) - (edge.xt - edge.xb);
 673
 674     // If new edge and ST edge don't cross
 675     if (edge.xt >= (*st).xt || edge.dx == (*st).dx || abs(den) <= DBL_EPSILON) {
 676       // No intersection - insert edge here (before the ST edge)
 677       existing_node = *st;
 678       *st = MALLOC!st_node();
 679       (*st).edge = edge;
 680       (*st).xb = edge.xb;
 681       (*st).xt = edge.xt;
 682       (*st).dx = edge.dx;
 683       (*st).prev = existing_node;
 684     } else {
 685       // Compute intersection between new edge and ST edge
 686       r = (edge.xb - (*st).xb) / den;
 687       x = (*st).xb + r * ((*st).xt - (*st).xb);
 688       y = r * dy;
 689
 690       // Insert the edge pointers and the intersection point in the IT
 691       add_intersection(it, (*st).edge, edge, x, y);
 692
 693       // Head further into the ST
 694       add_st_edge(&((*st).prev), it, edge, dy);
 695     }
 696   }
 697 }
 698
 699
 700 private void build_intersection_table (it_node** it, edge_node* aet, double dy) {
 701   st_node* st, stp;
 702   edge_node* edge;
 703
 704   // Build intersection table for the current scanbeam
 705   reset_it(it);
 706   st = null;
 707
 708   // Process each AET edge
 709   for (edge = aet; edge !is null; edge = edge.next) {
 710     if (edge.bstate[ABOVE] == BUNDLE_HEAD || edge.bundle[ABOVE][CLIP] || edge.bundle[ABOVE][SUBJ]) {
 711       add_st_edge(&st, it, edge, dy);
 712     }
 713   }
 714
 715   // Free the sorted edge table
 716   while (st !is null) {
 717     stp = st.prev;
 718     FREE(st);
 719     st = stp;
 720   }
 721 }
 722
 723
 724 private void swap_intersecting_edge_bundles (edge_node** aet, it_node* intersect) {
 725   edge_node* e0 = intersect.ie[0];
 726   edge_node* e1 = intersect.ie[1];
 727   edge_node* e0t = e0;
 728   edge_node* e1t = e1;
 729   edge_node* e0n = e0.next;
 730   edge_node* e1n = e1.next;
 731
 732   // Find the node before the e0 bundle
 733   edge_node* e0p = e0.prev;
 734   if (e0.bstate[ABOVE] == BUNDLE_HEAD) {
 735     do {
 736       e0t = e0p;
 737       e0p = e0p.prev;
 738     } while (e0p && (e0p.bstate[ABOVE] == BUNDLE_TAIL));
 739   }
 740
 741   // Find the node before the e1 bundle
 742   edge_node* e1p = e1.prev;
 743   if (e1.bstate[ABOVE] == BUNDLE_HEAD) {
 744     do {
 745       e1t = e1p;
 746       e1p = e1p.prev;
 747     } while (e1p && (e1p.bstate[ABOVE] == BUNDLE_TAIL));
 748   }
 749
 750   // Swap the e0p and e1p links
 751   if (e0p) {
 752     if (e1p) {
 753       if (e0p !is e1) {
 754         e0p.next = e1t;
 755         e1t.prev = e0p;
 756       }
 757       if (e1p !is e0) {
 758         e1p.next = e0t;
 759         e0t.prev = e1p;
 760       }
 761     } else {
 762       if (e0p !is e1) {
 763         e0p.next = e1t;
 764         e1t.prev = e0p;
 765       }
 766       *aet = e0t;
 767       e0t.prev = null;
 768     }
 769   } else {
 770     if (e1p !is e0) {
 771       e1p.next = e0t;
 772       e0t.prev = e1p;
 773     }
 774     *aet = e1t;
 775     e1t.prev = null;
 776   }
 777
 778   // Re-link after e0
 779   if (e0p !is e1) {
 780     e0.next = e1n;
 781     if (e1n) e1n.prev = e0;
 782   } else {
 783     e0.next = e1t;
 784     e1t.prev = e0;
 785   }
 786
 787   // Re-link after e1
 788   if (e1p !is e0) {
 789     e1.next = e0n;
 790     if (e0n) e0n.prev = e1;
 791   } else {
 792     e1.next = e0t;
 793     e0t.prev = e1;
 794   }
 795 }
 796
 797
 798 private int count_contours (polygon_node* polygon) {
 799   int nc, nv;
 800   vertex_node* v, nextv;
 801
 802   for (nc = 0; polygon !is null; polygon = polygon.next) {
 803     if (polygon.active) {
 804       // Count the vertices in the current contour
 805       nv = 0;
 806       for (v= polygon.proxy.v[LEFT]; v; v= v.next) ++nv;
 807
 808       // Record valid vertex counts in the active field
 809       if (nv > 2) {
 810         polygon.active= nv;
 811         ++nc;
 812       } else {
 813         // Invalid contour: just free the heap
 814         for (v = polygon.proxy.v[LEFT]; v !is null; v = nextv) {
 815           nextv = v.next;
 816           FREE(v);
 817         }
 818         polygon.active = 0;
 819       }
 820     }
 821   }
 822
 823   return nc;
 824 }
 825
 826
 827 private void add_left (polygon_node* p, double x, double y) {
 828   vertex_node* nv;
 829
 830   // Create a new vertex node and set its fields
 831   nv = MALLOC!vertex_node();
 832   nv.x = x;
 833   nv.y = y;
 834
 835   // Add vertex nv to the left end of the polygon's vertex list
 836   nv.next = p.proxy.v[LEFT];
 837
 838   // Update proxy.[LEFT] to point to nv
 839   p.proxy.v[LEFT] = nv;
 840 }
 841
 842
 843 private void merge_left (polygon_node* p, polygon_node* q, polygon_node* list) {
 844   polygon_node* target;
 845
 846   // Label contour as a hole
 847   q.proxy.hole = true;
 848
 849   if (p.proxy !is q.proxy) {
 850     // Assign p's vertex list to the left end of q's list
 851     p.proxy.v[RIGHT].next = q.proxy.v[LEFT];
 852     q.proxy.v[LEFT] = p.proxy.v[LEFT];
 853
 854     // Redirect any p.proxy references to q.proxy
 855     for (target = p.proxy; list !is null; list = list.next) {
 856       if (list.proxy is target) {
 857         list.active = false;
 858         list.proxy = q.proxy;
 859       }
 860     }
 861   }
 862 }
 863
 864
 865 private void add_right (polygon_node* p, double x, double y) {
 866   vertex_node* nv;
 867
 868   // Create a new vertex node and set its fields
 869   nv = MALLOC!vertex_node();
 870   nv.x = x;
 871   nv.y = y;
 872   nv.next = null;
 873
 874   // Add vertex nv to the right end of the polygon's vertex list
 875   p.proxy.v[RIGHT].next = nv;
 876
 877   // Update proxy.v[RIGHT] to point to nv
 878   p.proxy.v[RIGHT] = nv;
 879 }
 880
 881
 882 private void merge_right (polygon_node* p, polygon_node* q, polygon_node* list) {
 883   polygon_node* target;
 884
 885   // Label contour as external
 886   q.proxy.hole = false;
 887
 888   if (p.proxy !is q.proxy) {
 889     // Assign p's vertex list to the right end of q's list
 890     q.proxy.v[RIGHT].next = p.proxy.v[LEFT];
 891     q.proxy.v[RIGHT] = p.proxy.v[RIGHT];
 892
 893     // Redirect any p.proxy references to q.proxy
 894     for (target = p.proxy; list !is null; list = list.next) {
 895       if (list.proxy is target) {
 896         list.active = false;
 897         list.proxy = q.proxy;
 898       }
 899     }
 900   }
 901 }
 902
 903
 904 private void add_local_min (polygon_node** p, edge_node* edge, double x, double y) {
 905   polygon_node* existing_min;
 906   vertex_node* nv;
 907
 908   existing_min = *p;
 909
 910   *p = MALLOC!polygon_node();
 911
 912   // Create a new vertex node and set its fields
 913   nv = MALLOC!vertex_node();
 914   nv.x = x;
 915   nv.y = y;
 916   nv.next = null;
 917
 918   // Initialise proxy to point to p itself
 919   (*p).proxy = (*p);
 920   (*p).active = true;
 921   (*p).next = existing_min;
 922
 923   // Make v[LEFT] and v[RIGHT] point to new vertex nv
 924   (*p).v[LEFT] = nv;
 925   (*p).v[RIGHT] = nv;
 926
 927   // Assign polygon p to the edge
 928   edge.outp[ABOVE] = *p;
 929 }
 930
 931
 932 private int count_tristrips (polygon_node* tn) {
 933   int total;
 934   for (total = 0; tn !is null; tn = tn.next) if (tn.active > 2) ++total;
 935   return total;
 936 }
 937
 938
 939 private void add_vertex (vertex_node** t, double x, double y) {
 940   if (!(*t)) {
 941     *t = MALLOC!vertex_node();
 942     (*t).x = x;
 943     (*t).y = y;
 944     (*t).next = null;
 945   } else {
 946     // Head further down the list
 947     add_vertex(&((*t).next), x, y);
 948   }
 949 }
 950
 951
 952 private void new_tristrip (polygon_node** tn, edge_node* edge, double x, double y) {
 953   if (!(*tn)) {
 954     *tn = MALLOC!polygon_node();
 955     (*tn).next = null;
 956     (*tn).v[LEFT] = null;
 957     (*tn).v[RIGHT] = null;
 958     (*tn).active = 1;
 959     add_vertex(&((*tn).v[LEFT]), x, y);
 960     edge.outp[ABOVE] = *tn;
 961   } else {
 962     // Head further down the list
 963     new_tristrip(&((*tn).next), edge, x, y);
 964   }
 965 }
 966
 967
 968 private bbox* create_contour_bboxes (gpc_polygon* p) {
 969   bbox* box = MALLOC!bbox(p.num_contours);
 970
 971   // Construct contour bounding boxes
 972   foreach (immutable int c; 0..p.num_contours) {
 973     // Initialise bounding box extent
 974     box[c].xmin = double.max;
 975     box[c].ymin = double.max;
 976     box[c].xmax = -double.max;
 977     box[c].ymax = -double.max;
 978
 979     foreach (immutable int v; 0..p.contour[c].num_vertices) {
 980       // Adjust bounding box
 981       if (p.contour[c].vertex[v].x < box[c].xmin) box[c].xmin = p.contour[c].vertex[v].x;
 982       if (p.contour[c].vertex[v].y < box[c].ymin) box[c].ymin = p.contour[c].vertex[v].y;
 983       if (p.contour[c].vertex[v].x > box[c].xmax) box[c].xmax = p.contour[c].vertex[v].x;
 984       if (p.contour[c].vertex[v].y > box[c].ymax) box[c].ymax = p.contour[c].vertex[v].y;
 985     }
 986   }
 987
 988   return box;
 989 }
 990
 991
 992 private void minimax_test (gpc_polygon* subj, gpc_polygon* clip, GPC op) {
 993   bbox* s_bbox, c_bbox;
 994   int overlap;
 995   int* o_table;
 996
 997   s_bbox = create_contour_bboxes(subj);
 998   c_bbox = create_contour_bboxes(clip);
 999
1000   o_table = MALLOC!int(subj.num_contours*clip.num_contours);
1001
1002   // Check all subject contour bounding boxes against clip boxes
1003   foreach (immutable int s; 0..subj.num_contours) {
1004     foreach (immutable int c; 0..clip.num_contours) {
1005       o_table[c*subj.num_contours + s] =
1006              (!((s_bbox[s].xmax < c_bbox[c].xmin) ||
1007                 (s_bbox[s].xmin > c_bbox[c].xmax))) &&
1008              (!((s_bbox[s].ymax < c_bbox[c].ymin) ||
1009                 (s_bbox[s].ymin > c_bbox[c].ymax)));
1010     }
1011   }
1012
1013   // For each clip contour, search for any subject contour overlaps
1014   foreach (immutable int c; 0..clip.num_contours) {
1015     overlap = 0;
1016     for (int s = 0; !overlap && s < subj.num_contours; ++s) overlap = o_table[c * subj.num_contours + s];
1017     // Flag non contributing status by negating vertex count
1018     if (!overlap) clip.contour[c].num_vertices = -clip.contour[c].num_vertices;
1019   }
1020
1021   if (op == GPC.Int) {
1022     // For each subject contour, search for any clip contour overlaps
1023     foreach (immutable int s; 0..subj.num_contours) {
1024       overlap = 0;
1025       for (int c = 0; !overlap && c < clip.num_contours; ++c) overlap = o_table[c * subj.num_contours + s];
1026       // Flag non contributing status by negating vertex count
1027       if (!overlap) subj.contour[s].num_vertices = -subj.contour[s].num_vertices;
1028     }
1029   }
1030
1031   FREE(s_bbox);
1032   FREE(c_bbox);
1033   FREE(o_table);
1034 }
1035
1036
1037 /*
1038 ===========================================================================
1039                              Public Functions
1040 ===========================================================================
1041 */
1042
1043 import core.stdc.stdio : FILE;
1044
1045 public void gpc_read_polygon() (FILE* fp, bool read_hole_flags, ref gpc_polygon p) {
1046   import core.stdc.stdio : fscanf;
1047   fscanf(fp, "%d", &(p.num_contours));
1048   p.hole = MALLOC!int(p.num_contours);
1049   p.contour = MALLOC!gpc_vertex_list(p.num_contours);
1050   foreach (immutable int c; 0..p.num_contours) {
1051     fscanf(fp, "%d", &(p.contour[c].num_vertices));
1052     if (read_hole_flags) {
1053       fscanf(fp, "%d", &(p.hole[c]));
1054     } else {
1055       p.hole[c] = false; // Assume all contours to be external
1056     }
1057     p.contour[c].vertex = MALLOC!gpc_vertex(p.contour[c].num_vertices);
1058     foreach (immutable int v; 0..p.contour[c].num_vertices) {
1059       fscanf(fp, "%lf %lf", &(p.contour[c].vertex[v].x), &(p.contour[c].vertex[v].y));
1060     }
1061   }
1062 }
1063
1064 public void gpc_write_polygon() (FILE* fp, bool write_hole_flags, in ref gpc_polygon p) {
1065   import core.stdc.stdio : fprintf;
1066   fprintf(fp, "%d\n", p.num_contours);
1067   foreach (immutable int c; 0..p.num_contours) {
1068     fprintf(fp, "%d\n", p.contour[c].num_vertices);
1069     if (write_hole_flags) fprintf(fp, "%d\n", p.hole[c]);
1070     foreach (immutable int v; 0..p.contour[c].num_vertices) {
1071       fprintf(fp, "% .*lf % .*lf\n", double.dig, p.contour[c].vertex[v].x, double.dig, p.contour[c].vertex[v].y);
1072     }
1073   }
1074 }
1075
1076
1077 /// Frees allocated polygon memory. It is safe to call it on empty (but initialized) polygon.
1078 public void gpc_free_polygon (ref gpc_polygon p) {
1079   foreach (immutable int c; 0..p.num_contours) FREE(p.contour[c].vertex);
1080   FREE(p.hole);
1081   FREE(p.contour);
1082   p.num_contours = 0;
1083 }
1084
1085 /** Add new contour to the existing polygon.
1086  *
1087  * This function copies `new_contour` contents. But it reallocates quite inefficiently,
1088  * so you'd better build polys yourself.
1089  *
1090  * Contours with less than three points will not be added.
1091  *
1092  * Returns `true` if contour contains three or more points (i.e. it was added).
1093  */
1094 public bool gpc_add_contour (ref gpc_polygon p, const(gpc_vertex)[] new_contour, bool hole=false) {
1095   import core.stdc.string : memcpy;
1096
1097   if (new_contour.length < 3) return false;
1098   if (new_contour.length > int.max/8) assert(0, "out of memory");
1099
1100   int c = p.num_contours;
1101
1102   // allocate hole array with one more item
1103   REALLOC(p.hole, c+1);
1104
1105   // allocate contour array with one more item
1106   REALLOC(p.contour, c+1);
1107
1108   // copy the new contour and hole onto the end of the extended arrays
1109   p.hole[c] = hole;
1110   p.contour[c].num_vertices = cast(int)new_contour.length;
1111   p.contour[c].vertex = MALLOC!gpc_vertex(cast(uint)new_contour.length);
1112   memcpy(p.contour[c].vertex, new_contour.ptr, cast(uint)new_contour.length*gpc_vertex.sizeof);
1113
1114   // Update the polygon information
1115   ++p.num_contours;
1116
1117   return true;
1118 }
1119
1120
1121 /** Calculates clipping.
1122  *
1123  * [result] will be overwritten (so old contents won't be freed, and it may be uninitialized at all).
1124  */
1125 public void gpc_polygon_clip (GPC op, ref gpc_polygon subj, ref gpc_polygon clip, ref gpc_polygon result) {
1126   sb_tree* sbtree = null;
1127   it_node* it = null, intersect;
1128   edge_node* edge, prev_edge, next_edge, succ_edge, e0, e1;
1129   edge_node* aet = null, c_heap = null, s_heap = null;
1130   lmt_node* lmt = null, local_min;
1131   polygon_node* out_poly = null, p, q, poly, npoly, cf = null;
1132   vertex_node* vtx, nv;
1133   h_state[2] horiz;
1134   int[2] inn, exists;
1135   int[2] parity = [LEFT, LEFT];
1136   int c, v, contributing, scanbeam = 0, sbt_entries = 0;
1137   int vclass, bl, br, tl, tr;
1138   double* sbt = null;
1139   double xb, px, yb, yt, dy, ix, iy;
1140
1141   // Test for trivial NULL result cases
1142   if ((subj.num_contours == 0 && clip.num_contours == 0) ||
1143       (subj.num_contours == 0 && (op == GPC.Int || op == GPC.Diff)) ||
1144       (clip.num_contours == 0 && op == GPC.Int))
1145   {
1146     result.num_contours = 0;
1147     result.hole = null;
1148     result.contour = null;
1149     return;
1150   }
1151
1152   // Identify potentialy contributing contours
1153   if ((op == GPC.Int || op == GPC.Diff) && subj.num_contours > 0 && clip.num_contours > 0) minimax_test(&subj, &clip, op);
1154
1155   // Build LMT
1156   if (subj.num_contours > 0) s_heap = build_lmt(&lmt, &sbtree, &sbt_entries, &subj, SUBJ, op);
1157   if (clip.num_contours > 0) c_heap = build_lmt(&lmt, &sbtree, &sbt_entries, &clip, CLIP, op);
1158
1159   // Return a NULL result if no contours contribute
1160   if (lmt is null) {
1161     result.num_contours = 0;
1162     result.hole = null;
1163     result.contour = null;
1164     reset_lmt(&lmt);
1165     FREE(s_heap);
1166     FREE(c_heap);
1167     return;
1168   }
1169
1170   // Build scanbeam table from scanbeam tree
1171   sbt = MALLOC!double(sbt_entries);
1172   build_sbt(&scanbeam, sbt, sbtree);
1173   scanbeam = 0;
1174   free_sbtree(&sbtree);
1175
1176   // Allow pointer re-use without causing memory leak
1177   if (&subj is &result) gpc_free_polygon(subj);
1178   if (&clip is &result) gpc_free_polygon(clip);
1179
1180   // Invert clip polygon for difference operation
1181   if (op == GPC.Diff) parity[CLIP] = RIGHT;
1182
1183   local_min = lmt;
1184
1185   // Process each scanbeam
1186   while (scanbeam < sbt_entries) {
1187     // Set yb and yt to the bottom and top of the scanbeam
1188     yb = sbt[scanbeam++];
1189     if (scanbeam < sbt_entries) {
1190       yt = sbt[scanbeam];
1191       dy = yt - yb;
1192     }
1193
1194     // === SCANBEAM BOUNDARY PROCESSING ================================
1195
1196     // If LMT node corresponding to yb exists
1197     if (local_min) {
1198       if (local_min.y == yb) {
1199         // Add edges starting at this local minimum to the AET
1200         for (edge = local_min.first_bound; edge !is null; edge = edge.next_bound) add_edge_to_aet(&aet, edge, null);
1201         local_min = local_min.next;
1202       }
1203     }
1204
1205     // Set dummy previous x value
1206     px = -double.max;
1207
1208     // Create bundles within AET
1209     e0 = aet;
1210     e1 = aet;
1211
1212     // Set up bundle fields of first edge
1213     aet.bundle[ABOVE][ aet.type] = (aet.top.y != yb);
1214     aet.bundle[ABOVE][!aet.type] = false;
1215     aet.bstate[ABOVE] = UNBUNDLED;
1216
1217     for (next_edge = aet.next; next_edge !is null; next_edge = next_edge.next) {
1218       // Set up bundle fields of next edge
1219       next_edge.bundle[ABOVE][ next_edge.type] = (next_edge.top.y != yb);
1220       next_edge.bundle[ABOVE][!next_edge.type] = false;
1221       next_edge.bstate[ABOVE] = UNBUNDLED;
1222
1223       // Bundle edges above the scanbeam boundary if they coincide
1224       if (next_edge.bundle[ABOVE][next_edge.type]) {
1225         if (EQ(e0.xb, next_edge.xb) && EQ(e0.dx, next_edge.dx) && e0.top.y != yb) {
1226           next_edge.bundle[ABOVE][ next_edge.type] ^= e0.bundle[ABOVE][ next_edge.type];
1227           next_edge.bundle[ABOVE][!next_edge.type] = e0.bundle[ABOVE][!next_edge.type];
1228           next_edge.bstate[ABOVE] = BUNDLE_HEAD;
1229           e0.bundle[ABOVE][CLIP] = false;
1230           e0.bundle[ABOVE][SUBJ] = false;
1231           e0.bstate[ABOVE] = BUNDLE_TAIL;
1232         }
1233         e0 = next_edge;
1234       }
1235     }
1236
1237     horiz[CLIP] = NH;
1238     horiz[SUBJ] = NH;
1239
1240     // Process each edge at this scanbeam boundary
1241     for (edge = aet; edge !is null; edge = edge.next) {
1242       exists[CLIP] = edge.bundle[ABOVE][CLIP] + (edge.bundle[BELOW][CLIP] << 1);
1243       exists[SUBJ] = edge.bundle[ABOVE][SUBJ] + (edge.bundle[BELOW][SUBJ] << 1);
1244
1245       if (exists[CLIP] || exists[SUBJ]) {
1246         // Set bundle side
1247         edge.bside[CLIP] = parity[CLIP];
1248         edge.bside[SUBJ] = parity[SUBJ];
1249
1250         // Determine contributing status and quadrant occupancies
1251         switch (op) {
1252           case GPC.Diff:
1253           case GPC.Int:
1254             contributing = (exists[CLIP] && (parity[SUBJ] || horiz[SUBJ])) || (exists[SUBJ] && (parity[CLIP] || horiz[CLIP])) || (exists[CLIP] && exists[SUBJ] && (parity[CLIP] == parity[SUBJ]));
1255             br = (parity[CLIP]) && (parity[SUBJ]);
1256             bl = (parity[CLIP] ^ edge.bundle[ABOVE][CLIP]) && (parity[SUBJ] ^ edge.bundle[ABOVE][SUBJ]);
1257             tr = (parity[CLIP] ^ (horiz[CLIP] != NH)) && (parity[SUBJ] ^ (horiz[SUBJ] != NH));
1258             tl = (parity[CLIP] ^ (horiz[CLIP] != NH) ^ edge.bundle[BELOW][CLIP]) && (parity[SUBJ] ^ (horiz[SUBJ] != NH) ^ edge.bundle[BELOW][SUBJ]);
1259             break;
1260           case GPC.Xor:
1261             contributing = exists[CLIP] || exists[SUBJ];
1262             br = (parity[CLIP]) ^ (parity[SUBJ]);
1263             bl = (parity[CLIP] ^ edge.bundle[ABOVE][CLIP]) ^ (parity[SUBJ] ^ edge.bundle[ABOVE][SUBJ]);
1264             tr = (parity[CLIP] ^ (horiz[CLIP] != NH)) ^ (parity[SUBJ] ^ (horiz[SUBJ] != NH));
1265             tl = (parity[CLIP] ^ (horiz[CLIP] != NH) ^ edge.bundle[BELOW][CLIP]) ^ (parity[SUBJ] ^ (horiz[SUBJ] != NH) ^ edge.bundle[BELOW][SUBJ]);
1266             break;
1267           case GPC.Union:
1268             contributing = (exists[CLIP] && (!parity[SUBJ] || horiz[SUBJ])) || (exists[SUBJ] && (!parity[CLIP] || horiz[CLIP])) || (exists[CLIP] && exists[SUBJ] && (parity[CLIP] == parity[SUBJ]));
1269             br = (parity[CLIP]) || (parity[SUBJ]);
1270             bl = (parity[CLIP] ^ edge.bundle[ABOVE][CLIP]) || (parity[SUBJ] ^ edge.bundle[ABOVE][SUBJ]);
1271             tr = (parity[CLIP] ^ (horiz[CLIP] != NH)) || (parity[SUBJ] ^ (horiz[SUBJ] != NH));
1272             tl = (parity[CLIP] ^ (horiz[CLIP] != NH) ^ edge.bundle[BELOW][CLIP]) || (parity[SUBJ] ^ (horiz[SUBJ] != NH) ^ edge.bundle[BELOW][SUBJ]);
1273             break;
1274           default: assert(0);
1275         }
1276
1277         // Update parity
1278         parity[CLIP] ^= edge.bundle[ABOVE][CLIP];
1279         parity[SUBJ] ^= edge.bundle[ABOVE][SUBJ];
1280
1281         // Update horizontal state
1282         if (exists[CLIP]) horiz[CLIP] = next_h_state[horiz[CLIP]][((exists[CLIP] - 1) << 1) + parity[CLIP]];
1283         if (exists[SUBJ]) horiz[SUBJ] = next_h_state[horiz[SUBJ]][((exists[SUBJ] - 1) << 1) + parity[SUBJ]];
1284
1285         vclass = tr + (tl << 1) + (br << 2) + (bl << 3);
1286
1287         if (contributing) {
1288           xb = edge.xb;
1289
1290           switch (vclass) {
1291             case EMN:
1292             case IMN:
1293               add_local_min(&out_poly, edge, xb, yb);
1294               px = xb;
1295               cf = edge.outp[ABOVE];
1296               break;
1297             case ERI:
1298               if (xb != px) {
1299                 add_right(cf, xb, yb);
1300                 px = xb;
1301               }
1302               edge.outp[ABOVE]= cf;
1303               cf = null;
1304               break;
1305             case ELI:
1306               add_left(edge.outp[BELOW], xb, yb);
1307               px = xb;
1308               cf = edge.outp[BELOW];
1309               break;
1310             case EMX:
1311               if (xb != px) {
1312                 add_left(cf, xb, yb);
1313                 px = xb;
1314               }
1315               merge_right(cf, edge.outp[BELOW], out_poly);
1316               cf = null;
1317               break;
1318             case ILI:
1319               if (xb != px) {
1320                 add_left(cf, xb, yb);
1321                 px = xb;
1322               }
1323               edge.outp[ABOVE] = cf;
1324               cf = null;
1325               break;
1326             case IRI:
1327               add_right(edge.outp[BELOW], xb, yb);
1328               px = xb;
1329               cf = edge.outp[BELOW];
1330               edge.outp[BELOW] = null;
1331               break;
1332             case IMX:
1333               if (xb != px) {
1334                 add_right(cf, xb, yb);
1335                 px = xb;
1336               }
1337               merge_left(cf, edge.outp[BELOW], out_poly);
1338               cf = null;
1339               edge.outp[BELOW] = null;
1340               break;
1341             case IMM:
1342               if (xb != px) {
1343                 add_right(cf, xb, yb);
1344                 px = xb;
1345               }
1346               merge_left(cf, edge.outp[BELOW], out_poly);
1347               edge.outp[BELOW] = null;
1348               add_local_min(&out_poly, edge, xb, yb);
1349               cf = edge.outp[ABOVE];
1350               break;
1351             case EMM:
1352               if (xb != px) {
1353                 add_left(cf, xb, yb);
1354                 px = xb;
1355               }
1356               merge_right(cf, edge.outp[BELOW], out_poly);
1357               edge.outp[BELOW] = null;
1358               add_local_min(&out_poly, edge, xb, yb);
1359               cf = edge.outp[ABOVE];
1360               break;
1361             case LED:
1362               if (edge.bot.y == yb) add_left(edge.outp[BELOW], xb, yb);
1363               edge.outp[ABOVE] = edge.outp[BELOW];
1364               px = xb;
1365               break;
1366             case RED:
1367               if (edge.bot.y == yb) add_right(edge.outp[BELOW], xb, yb);
1368               edge.outp[ABOVE] = edge.outp[BELOW];
1369               px = xb;
1370               break;
1371             default:
1372               break;
1373           } // End of switch
1374         } // End of contributing conditional
1375       } // End of edge exists conditional
1376     } // End of AET loop
1377
1378     // Delete terminating edges from the AET, otherwise compute xt
1379     for (edge = aet; edge !is null; edge = edge.next) {
1380       if (edge.top.y == yb) {
1381         prev_edge = edge.prev;
1382         next_edge = edge.next;
1383         if (prev_edge) prev_edge.next = next_edge; else aet = next_edge;
1384         if (next_edge) next_edge.prev = prev_edge;
1385
1386         // Copy bundle head state to the adjacent tail edge if required
1387         if ((edge.bstate[BELOW] == BUNDLE_HEAD) && prev_edge) {
1388           if (prev_edge.bstate[BELOW] == BUNDLE_TAIL) {
1389             prev_edge.outp[BELOW] = edge.outp[BELOW];
1390             prev_edge.bstate[BELOW] = UNBUNDLED;
1391             if (prev_edge.prev) {
1392               if (prev_edge.prev.bstate[BELOW] == BUNDLE_TAIL) prev_edge.bstate[BELOW] = BUNDLE_HEAD;
1393             }
1394           }
1395         }
1396       } else {
1397         if (edge.top.y == yt) {
1398           edge.xt = edge.top.x;
1399         } else {
1400           edge.xt = edge.bot.x + edge.dx * (yt - edge.bot.y);
1401         }
1402       }
1403     }
1404
1405     if (scanbeam < sbt_entries) {
1406       // === SCANBEAM INTERIOR PROCESSING ==============================
1407
1408       build_intersection_table(&it, aet, dy);
1409
1410       // Process each node in the intersection table
1411       for (intersect = it; intersect; intersect = intersect.next) {
1412         e0 = intersect.ie[0];
1413         e1 = intersect.ie[1];
1414
1415         // Only generate output for contributing intersections
1416         if ((e0.bundle[ABOVE][CLIP] || e0.bundle[ABOVE][SUBJ]) && (e1.bundle[ABOVE][CLIP] || e1.bundle[ABOVE][SUBJ])) {
1417           p = e0.outp[ABOVE];
1418           q = e1.outp[ABOVE];
1419           ix = intersect.point.x;
1420           iy = intersect.point.y + yb;
1421
1422           inn[CLIP] = (e0.bundle[ABOVE][CLIP] && !e0.bside[CLIP]) || (e1.bundle[ABOVE][CLIP] && e1.bside[CLIP]) || (!e0.bundle[ABOVE][CLIP] && !e1.bundle[ABOVE][CLIP] && e0.bside[CLIP] && e1.bside[CLIP]);
1423           inn[SUBJ] = (e0.bundle[ABOVE][SUBJ] && !e0.bside[SUBJ]) || (e1.bundle[ABOVE][SUBJ] && e1.bside[SUBJ]) || (!e0.bundle[ABOVE][SUBJ] && !e1.bundle[ABOVE][SUBJ] && e0.bside[SUBJ] && e1.bside[SUBJ]);
1424
1425           // Determine quadrant occupancies
1426           switch (op) {
1427             case GPC.Diff:
1428             case GPC.Int:
1429               tr = (inn[CLIP]) && (inn[SUBJ]);
1430               tl = (inn[CLIP] ^ e1.bundle[ABOVE][CLIP]) && (inn[SUBJ] ^ e1.bundle[ABOVE][SUBJ]);
1431               br = (inn[CLIP] ^ e0.bundle[ABOVE][CLIP]) && (inn[SUBJ] ^ e0.bundle[ABOVE][SUBJ]);
1432               bl = (inn[CLIP] ^ e1.bundle[ABOVE][CLIP] ^ e0.bundle[ABOVE][CLIP]) && (inn[SUBJ] ^ e1.bundle[ABOVE][SUBJ] ^ e0.bundle[ABOVE][SUBJ]);
1433               break;
1434             case GPC.Xor:
1435               tr = (inn[CLIP]) ^ (inn[SUBJ]);
1436               tl = (inn[CLIP] ^ e1.bundle[ABOVE][CLIP]) ^ (inn[SUBJ] ^ e1.bundle[ABOVE][SUBJ]);
1437               br = (inn[CLIP] ^ e0.bundle[ABOVE][CLIP]) ^ (inn[SUBJ] ^ e0.bundle[ABOVE][SUBJ]);
1438               bl = (inn[CLIP] ^ e1.bundle[ABOVE][CLIP] ^ e0.bundle[ABOVE][CLIP]) ^ (inn[SUBJ] ^ e1.bundle[ABOVE][SUBJ] ^ e0.bundle[ABOVE][SUBJ]);
1439               break;
1440             case GPC.Union:
1441               tr = (inn[CLIP]) || (inn[SUBJ]);
1442               tl = (inn[CLIP] ^ e1.bundle[ABOVE][CLIP]) || (inn[SUBJ] ^ e1.bundle[ABOVE][SUBJ]);
1443               br = (inn[CLIP] ^ e0.bundle[ABOVE][CLIP]) || (inn[SUBJ] ^ e0.bundle[ABOVE][SUBJ]);
1444               bl = (inn[CLIP] ^ e1.bundle[ABOVE][CLIP] ^ e0.bundle[ABOVE][CLIP]) || (inn[SUBJ] ^ e1.bundle[ABOVE][SUBJ] ^ e0.bundle[ABOVE][SUBJ]);
1445               break;
1446             default: assert(0);
1447           }
1448
1449           vclass = tr + (tl << 1) + (br << 2) + (bl << 3);
1450
1451           switch (vclass) {
1452             case EMN:
1453               add_local_min(&out_poly, e0, ix, iy);
1454               e1.outp[ABOVE] = e0.outp[ABOVE];
1455               break;
1456             case ERI:
1457               if (p) {
1458                 add_right(p, ix, iy);
1459                 e1.outp[ABOVE] = p;
1460                 e0.outp[ABOVE] = null;
1461               }
1462               break;
1463             case ELI:
1464               if (q) {
1465                 add_left(q, ix, iy);
1466                 e0.outp[ABOVE] = q;
1467                 e1.outp[ABOVE] = null;
1468               }
1469               break;
1470             case EMX:
1471               if (p && q) {
1472                 add_left(p, ix, iy);
1473                 merge_right(p, q, out_poly);
1474                 e0.outp[ABOVE] = null;
1475                 e1.outp[ABOVE] = null;
1476               }
1477               break;
1478             case IMN:
1479               add_local_min(&out_poly, e0, ix, iy);
1480               e1.outp[ABOVE] = e0.outp[ABOVE];
1481               break;
1482             case ILI:
1483               if (p) {
1484                 add_left(p, ix, iy);
1485                 e1.outp[ABOVE] = p;
1486                 e0.outp[ABOVE] = null;
1487               }
1488               break;
1489             case IRI:
1490               if (q) {
1491                 add_right(q, ix, iy);
1492                 e0.outp[ABOVE] = q;
1493                 e1.outp[ABOVE] = null;
1494               }
1495               break;
1496             case IMX:
1497               if (p && q) {
1498                 add_right(p, ix, iy);
1499                 merge_left(p, q, out_poly);
1500                 e0.outp[ABOVE] = null;
1501                 e1.outp[ABOVE] = null;
1502               }
1503               break;
1504             case IMM:
1505               if (p && q) {
1506                 add_right(p, ix, iy);
1507                 merge_left(p, q, out_poly);
1508                 add_local_min(&out_poly, e0, ix, iy);
1509                 e1.outp[ABOVE] = e0.outp[ABOVE];
1510               }
1511               break;
1512             case EMM:
1513               if (p && q) {
1514                 add_left(p, ix, iy);
1515                 merge_right(p, q, out_poly);
1516                 add_local_min(&out_poly, e0, ix, iy);
1517                 e1.outp[ABOVE] = e0.outp[ABOVE];
1518               }
1519               break;
1520             default:
1521               break;
1522           } // End of switch
1523         } // End of contributing intersection conditional
1524
1525         // Swap bundle sides in response to edge crossing
1526         if (e0.bundle[ABOVE][CLIP]) e1.bside[CLIP] = !e1.bside[CLIP];
1527         if (e1.bundle[ABOVE][CLIP]) e0.bside[CLIP] = !e0.bside[CLIP];
1528         if (e0.bundle[ABOVE][SUBJ]) e1.bside[SUBJ] = !e1.bside[SUBJ];
1529         if (e1.bundle[ABOVE][SUBJ]) e0.bside[SUBJ] = !e0.bside[SUBJ];
1530
1531         // Swap the edge bundles in the aet
1532         swap_intersecting_edge_bundles(&aet, intersect);
1533       } // End of IT loop
1534
1535       // Prepare for next scanbeam
1536       for (edge = aet; edge !is null; edge = next_edge) {
1537         next_edge = edge.next;
1538         succ_edge = edge.succ;
1539
1540         if (edge.top.y == yt && succ_edge) {
1541           // Replace AET edge by its successor
1542           succ_edge.outp[BELOW] = edge.outp[ABOVE];
1543           succ_edge.bstate[BELOW] = edge.bstate[ABOVE];
1544           succ_edge.bundle[BELOW][CLIP] = edge.bundle[ABOVE][CLIP];
1545           succ_edge.bundle[BELOW][SUBJ] = edge.bundle[ABOVE][SUBJ];
1546           prev_edge = edge.prev;
1547           if (prev_edge) prev_edge.next = succ_edge; else aet = succ_edge;
1548           if (next_edge) next_edge.prev= succ_edge;
1549           succ_edge.prev = prev_edge;
1550           succ_edge.next = next_edge;
1551         } else {
1552           // Update this edge
1553           edge.outp[BELOW] = edge.outp[ABOVE];
1554           edge.bstate[BELOW] = edge.bstate[ABOVE];
1555           edge.bundle[BELOW][CLIP] = edge.bundle[ABOVE][CLIP];
1556           edge.bundle[BELOW][SUBJ] = edge.bundle[ABOVE][SUBJ];
1557           edge.xb = edge.xt;
1558         }
1559         edge.outp[ABOVE] = null;
1560       }
1561     }
1562   } // === END OF SCANBEAM PROCESSING ==================================
1563
1564   // Generate result polygon from out_poly
1565   result.contour = null;
1566   result.hole = null;
1567   result.num_contours = count_contours(out_poly);
1568   if (result.num_contours > 0) {
1569     result.hole = MALLOC!int(result.num_contours);
1570     result.contour = MALLOC!gpc_vertex_list(result.num_contours);
1571     c = 0;
1572     for (poly = out_poly; poly !is null; poly = npoly) {
1573       npoly = poly.next;
1574       if (poly.active) {
1575         result.hole[c] = poly.proxy.hole;
1576         result.contour[c].num_vertices = poly.active;
1577         result.contour[c].vertex = MALLOC!gpc_vertex(result.contour[c].num_vertices);
1578         v = result.contour[c].num_vertices-1;
1579         for (vtx = poly.proxy.v[LEFT]; vtx !is null; vtx = nv) {
1580           nv = vtx.next;
1581           result.contour[c].vertex[v].x = vtx.x;
1582           result.contour[c].vertex[v].y = vtx.y;
1583           FREE(vtx);
1584           --v;
1585         }
1586         ++c;
1587       }
1588       FREE(poly);
1589     }
1590   } else {
1591     for (poly = out_poly; poly !is null; poly = npoly) {
1592       npoly = poly.next;
1593       FREE(poly);
1594     }
1595   }
1596
1597   // Tidy up
1598   reset_it(&it);
1599   reset_lmt(&lmt);
1600   FREE(c_heap);
1601   FREE(s_heap);
1602   FREE(sbt);
1603 }
1604
1605
1606 /// Frees allocated tristrip memory. It is safe to call it on empty (but initialized) tristrip.
1607 public void gpc_free_tristrip (ref gpc_tristrip t) {
1608   foreach (immutable int s; 0..t.num_strips) FREE(t.strip[s].vertex);
1609   FREE(t.strip);
1610   t.num_strips = 0;
1611 }
1612
1613
1614 /** Converts polygon to triangle strip (tristrip).
1615  *
1616  * [result] will be overwritten (so old contents won't be freed, and it may be uninitialized at all).
1617  */
1618 public void gpc_polygon_to_tristrip (ref gpc_polygon s, ref gpc_tristrip result) {
1619   gpc_polygon c;
1620   c.num_contours = 0;
1621   c.hole = null;
1622   c.contour = null;
1623   gpc_tristrip_clip(GPC.Diff, s, c, result);
1624 }
1625
1626
1627 /** Calculates clipping, returns result as triangle strip (tristrip).
1628  *
1629  * [result] will be overwritten (so old contents won't be freed, and it may be uninitialized at all).
1630  */
1631 public void gpc_tristrip_clip (GPC op, ref gpc_polygon subj, ref gpc_polygon clip, ref gpc_tristrip result) {
1632   sb_tree* sbtree = null;
1633   it_node* it = null, intersect;
1634   edge_node* edge, prev_edge, next_edge, succ_edge, e0, e1;
1635   edge_node* aet = null, c_heap = null, s_heap = null, cf;
1636   lmt_node* lmt = null, local_min;
1637   polygon_node* tlist = null, tn, tnn, p, q;
1638   vertex_node* lt, ltn, rt, rtn;
1639   h_state[2] horiz;
1640   vertex_type cft;
1641   int[2] inn, exists;
1642   int[2] parity = [LEFT, LEFT];
1643   int s, v, contributing, scanbeam = 0, sbt_entries = 0;
1644   int vclass, bl, br, tl, tr;
1645   double* sbt = null;
1646   double xb, px, nx, yb, yt, dy, ix, iy;
1647
1648   // Test for trivial NULL result cases
1649   if ((subj.num_contours == 0 && clip.num_contours == 0) ||
1650       (subj.num_contours == 0 && (op == GPC.Int || op == GPC.Diff)) ||
1651       (clip.num_contours == 0 && op == GPC.Int))
1652   {
1653     result.num_strips = 0;
1654     result.strip = null;
1655     return;
1656   }
1657
1658   // Identify potentialy contributing contours
1659   if ((op == GPC.Int || op == GPC.Diff) && subj.num_contours > 0 && clip.num_contours > 0) minimax_test(&subj, &clip, op);
1660
1661   // Build LMT
1662   if (subj.num_contours > 0) s_heap = build_lmt(&lmt, &sbtree, &sbt_entries, &subj, SUBJ, op);
1663   if (clip.num_contours > 0) c_heap = build_lmt(&lmt, &sbtree, &sbt_entries, &clip, CLIP, op);
1664
1665   // Return a NULL result if no contours contribute
1666   if (lmt is null) {
1667     result.num_strips = 0;
1668     result.strip = null;
1669     reset_lmt(&lmt);
1670     FREE(s_heap);
1671     FREE(c_heap);
1672     return;
1673   }
1674
1675   // Build scanbeam table from scanbeam tree
1676   sbt = MALLOC!double(sbt_entries);
1677   build_sbt(&scanbeam, sbt, sbtree);
1678   scanbeam = 0;
1679   free_sbtree(&sbtree);
1680
1681   // Invert clip polygon for difference operation
1682   if (op == GPC.Diff) parity[CLIP] = RIGHT;
1683
1684   local_min = lmt;
1685
1686   // Process each scanbeam
1687   while (scanbeam < sbt_entries) {
1688     // Set yb and yt to the bottom and top of the scanbeam
1689     yb = sbt[scanbeam++];
1690     if (scanbeam < sbt_entries) {
1691       yt = sbt[scanbeam];
1692       dy = yt - yb;
1693     }
1694
1695     // === SCANBEAM BOUNDARY PROCESSING ================================
1696
1697     // If LMT node corresponding to yb exists
1698     if (local_min) {
1699       if (local_min.y == yb) {
1700         // Add edges starting at this local minimum to the AET
1701         for (edge = local_min.first_bound; edge !is null; edge = edge.next_bound) add_edge_to_aet(&aet, edge, null);
1702         local_min = local_min.next;
1703       }
1704     }
1705
1706     // Set dummy previous x value
1707     px = -double.max;
1708
1709     // Create bundles within AET
1710     e0 = aet;
1711     e1 = aet;
1712
1713     // Set up bundle fields of first edge
1714     aet.bundle[ABOVE][ aet.type] = (aet.top.y != yb);
1715     aet.bundle[ABOVE][!aet.type] = false;
1716     aet.bstate[ABOVE] = UNBUNDLED;
1717
1718     for (next_edge = aet.next; next_edge !is null; next_edge = next_edge.next) {
1719       // Set up bundle fields of next edge
1720       next_edge.bundle[ABOVE][next_edge.type] = (next_edge.top.y != yb);
1721       next_edge.bundle[ABOVE][!next_edge.type] = false;
1722       next_edge.bstate[ABOVE] = UNBUNDLED;
1723
1724       // Bundle edges above the scanbeam boundary if they coincide
1725       if (next_edge.bundle[ABOVE][next_edge.type]) {
1726         if (EQ(e0.xb, next_edge.xb) && EQ(e0.dx, next_edge.dx) && e0.top.y != yb) {
1727           next_edge.bundle[ABOVE][ next_edge.type] ^= e0.bundle[ABOVE][ next_edge.type];
1728           next_edge.bundle[ABOVE][!next_edge.type] = e0.bundle[ABOVE][!next_edge.type];
1729           next_edge.bstate[ABOVE] = BUNDLE_HEAD;
1730           e0.bundle[ABOVE][CLIP] = false;
1731           e0.bundle[ABOVE][SUBJ] = false;
1732           e0.bstate[ABOVE] = BUNDLE_TAIL;
1733         }
1734         e0 = next_edge;
1735       }
1736     }
1737
1738     horiz[CLIP] = NH;
1739     horiz[SUBJ] = NH;
1740
1741     // Process each edge at this scanbeam boundary
1742     for (edge = aet; edge !is null; edge = edge.next) {
1743       exists[CLIP] = edge.bundle[ABOVE][CLIP] + (edge.bundle[BELOW][CLIP] << 1);
1744       exists[SUBJ] = edge.bundle[ABOVE][SUBJ] + (edge.bundle[BELOW][SUBJ] << 1);
1745
1746       if (exists[CLIP] || exists[SUBJ]) {
1747         // Set bundle side
1748         edge.bside[CLIP] = parity[CLIP];
1749         edge.bside[SUBJ] = parity[SUBJ];
1750
1751         // Determine contributing status and quadrant occupancies
1752         switch (op) {
1753           case GPC.Diff:
1754           case GPC.Int:
1755             contributing= (exists[CLIP] && (parity[SUBJ] || horiz[SUBJ])) || (exists[SUBJ] && (parity[CLIP] || horiz[CLIP])) || (exists[CLIP] && exists[SUBJ] && (parity[CLIP] == parity[SUBJ]));
1756             br = (parity[CLIP]) && (parity[SUBJ]);
1757             bl = (parity[CLIP] ^ edge.bundle[ABOVE][CLIP]) && (parity[SUBJ] ^ edge.bundle[ABOVE][SUBJ]);
1758             tr = (parity[CLIP] ^ (horiz[CLIP] != NH)) && (parity[SUBJ] ^ (horiz[SUBJ] != NH));
1759             tl = (parity[CLIP] ^ (horiz[CLIP] != NH) ^ edge.bundle[BELOW][CLIP]) && (parity[SUBJ] ^ (horiz[SUBJ] != NH) ^ edge.bundle[BELOW][SUBJ]);
1760             break;
1761           case GPC.Xor:
1762             contributing = exists[CLIP] || exists[SUBJ];
1763             br = (parity[CLIP]) ^ (parity[SUBJ]);
1764             bl = (parity[CLIP] ^ edge.bundle[ABOVE][CLIP]) ^ (parity[SUBJ] ^ edge.bundle[ABOVE][SUBJ]);
1765             tr = (parity[CLIP] ^ (horiz[CLIP] != NH)) ^ (parity[SUBJ] ^ (horiz[SUBJ] != NH));
1766             tl = (parity[CLIP] ^ (horiz[CLIP] != NH) ^ edge.bundle[BELOW][CLIP]) ^ (parity[SUBJ] ^ (horiz[SUBJ] != NH) ^ edge.bundle[BELOW][SUBJ]);
1767             break;
1768           case GPC.Union:
1769             contributing = (exists[CLIP] && (!parity[SUBJ] || horiz[SUBJ])) || (exists[SUBJ] && (!parity[CLIP] || horiz[CLIP])) || (exists[CLIP] && exists[SUBJ] && (parity[CLIP] == parity[SUBJ]));
1770             br = (parity[CLIP]) || (parity[SUBJ]);
1771             bl = (parity[CLIP] ^ edge.bundle[ABOVE][CLIP]) || (parity[SUBJ] ^ edge.bundle[ABOVE][SUBJ]);
1772             tr = (parity[CLIP] ^ (horiz[CLIP] != NH)) || (parity[SUBJ] ^ (horiz[SUBJ] != NH));
1773             tl = (parity[CLIP] ^ (horiz[CLIP] != NH) ^ edge.bundle[BELOW][CLIP]) || (parity[SUBJ] ^ (horiz[SUBJ] != NH) ^ edge.bundle[BELOW][SUBJ]);
1774             break;
1775           default: assert(0);
1776         }
1777
1778         // Update parity
1779         parity[CLIP] ^= edge.bundle[ABOVE][CLIP];
1780         parity[SUBJ] ^= edge.bundle[ABOVE][SUBJ];
1781
1782         // Update horizontal state
1783         if (exists[CLIP]) horiz[CLIP] = next_h_state[horiz[CLIP]][((exists[CLIP] - 1) << 1) + parity[CLIP]];
1784         if (exists[SUBJ]) horiz[SUBJ] = next_h_state[horiz[SUBJ]][((exists[SUBJ] - 1) << 1) + parity[SUBJ]];
1785
1786         vclass = tr + (tl << 1) + (br << 2) + (bl << 3);
1787
1788         if (contributing) {
1789           xb = edge.xb;
1790
1791           switch (vclass) {
1792             case EMN:
1793               new_tristrip(&tlist, edge, xb, yb);
1794               cf = edge;
1795               break;
1796             case ERI:
1797               edge.outp[ABOVE] = cf.outp[ABOVE];
1798               if (xb != cf.xb) VERTEX(edge, ABOVE, RIGHT, xb, yb);
1799               cf = null;
1800               break;
1801             case ELI:
1802               VERTEX(edge, BELOW, LEFT, xb, yb);
1803               edge.outp[ABOVE] = null;
1804               cf = edge;
1805               break;
1806             case EMX:
1807               if (xb != cf.xb) VERTEX(edge, BELOW, RIGHT, xb, yb);
1808               edge.outp[ABOVE] = null;
1809               cf = null;
1810               break;
1811             case IMN:
1812               if (cft == LED) {
1813                 if (cf.bot.y != yb) VERTEX(cf, BELOW, LEFT, cf.xb, yb);
1814                 new_tristrip(&tlist, cf, cf.xb, yb);
1815               }
1816               edge.outp[ABOVE] = cf.outp[ABOVE];
1817               VERTEX(edge, ABOVE, RIGHT, xb, yb);
1818               break;
1819             case ILI:
1820               new_tristrip(&tlist, edge, xb, yb);
1821               cf = edge;
1822               cft = ILI;
1823               break;
1824             case IRI:
1825               if (cft == LED) {
1826                 if (cf.bot.y != yb) VERTEX(cf, BELOW, LEFT, cf.xb, yb);
1827                 new_tristrip(&tlist, cf, cf.xb, yb);
1828               }
1829               VERTEX(edge, BELOW, RIGHT, xb, yb);
1830               edge.outp[ABOVE] = null;
1831               break;
1832             case IMX:
1833               VERTEX(edge, BELOW, LEFT, xb, yb);
1834               edge.outp[ABOVE] = null;
1835               cft = IMX;
1836               break;
1837             case IMM:
1838               VERTEX(edge, BELOW, LEFT, xb, yb);
1839               edge.outp[ABOVE] = cf.outp[ABOVE];
1840               if (xb != cf.xb) VERTEX(cf, ABOVE, RIGHT, xb, yb);
1841               cf = edge;
1842               break;
1843             case EMM:
1844               VERTEX(edge, BELOW, RIGHT, xb, yb);
1845               edge.outp[ABOVE] = null;
1846               new_tristrip(&tlist, edge, xb, yb);
1847               cf = edge;
1848               break;
1849             case LED:
1850               if (edge.bot.y == yb) VERTEX(edge, BELOW, LEFT, xb, yb);
1851               edge.outp[ABOVE] = edge.outp[BELOW];
1852               cf = edge;
1853               cft = LED;
1854               break;
1855             case RED:
1856               edge.outp[ABOVE]= cf.outp[ABOVE];
1857               if (cft == LED) {
1858                 if (cf.bot.y == yb) {
1859                   VERTEX(edge, BELOW, RIGHT, xb, yb);
1860                 } else {
1861                   if (edge.bot.y == yb) {
1862                     VERTEX(cf, BELOW, LEFT, cf.xb, yb);
1863                     VERTEX(edge, BELOW, RIGHT, xb, yb);
1864                   }
1865                 }
1866               } else {
1867                 VERTEX(edge, BELOW, RIGHT, xb, yb);
1868                 VERTEX(edge, ABOVE, RIGHT, xb, yb);
1869               }
1870               cf = null;
1871               break;
1872             default:
1873               break;
1874           } // End of switch
1875         } // End of contributing conditional
1876       } // End of edge exists conditional
1877     } // End of AET loop
1878
1879     // Delete terminating edges from the AET, otherwise compute xt
1880     for (edge = aet; edge !is null; edge = edge.next) {
1881       if (edge.top.y == yb) {
1882         prev_edge = edge.prev;
1883         next_edge = edge.next;
1884         if (prev_edge) prev_edge.next = next_edge; else aet = next_edge;
1885         if (next_edge) next_edge.prev = prev_edge;
1886
1887         // Copy bundle head state to the adjacent tail edge if required
1888         if ((edge.bstate[BELOW] == BUNDLE_HEAD) && prev_edge) {
1889           if (prev_edge.bstate[BELOW] == BUNDLE_TAIL) {
1890             prev_edge.outp[BELOW] = edge.outp[BELOW];
1891             prev_edge.bstate[BELOW] = UNBUNDLED;
1892             if (prev_edge.prev) {
1893               if (prev_edge.prev.bstate[BELOW] == BUNDLE_TAIL) prev_edge.bstate[BELOW] = BUNDLE_HEAD;
1894             }
1895           }
1896         }
1897       } else {
1898         if (edge.top.y == yt) edge.xt = edge.top.x; else edge.xt = edge.bot.x + edge.dx * (yt - edge.bot.y);
1899       }
1900     }
1901
1902     if (scanbeam < sbt_entries) {
1903       // === SCANBEAM INTERIOR PROCESSING ==============================
1904
1905       build_intersection_table(&it, aet, dy);
1906
1907       // Process each node in the intersection table
1908       for (intersect = it; intersect !is null; intersect = intersect.next) {
1909         e0 = intersect.ie[0];
1910         e1 = intersect.ie[1];
1911
1912         // Only generate output for contributing intersections
1913         if ((e0.bundle[ABOVE][CLIP] || e0.bundle[ABOVE][SUBJ]) && (e1.bundle[ABOVE][CLIP] || e1.bundle[ABOVE][SUBJ])) {
1914           p = e0.outp[ABOVE];
1915           q = e1.outp[ABOVE];
1916           ix = intersect.point.x;
1917           iy = intersect.point.y + yb;
1918
1919           inn[CLIP] = (e0.bundle[ABOVE][CLIP] && !e0.bside[CLIP]) || (e1.bundle[ABOVE][CLIP] && e1.bside[CLIP]) || (!e0.bundle[ABOVE][CLIP] && !e1.bundle[ABOVE][CLIP] && e0.bside[CLIP] && e1.bside[CLIP]);
1920           inn[SUBJ] = (e0.bundle[ABOVE][SUBJ] && !e0.bside[SUBJ]) || (e1.bundle[ABOVE][SUBJ] && e1.bside[SUBJ]) || (!e0.bundle[ABOVE][SUBJ] && !e1.bundle[ABOVE][SUBJ] && e0.bside[SUBJ] && e1.bside[SUBJ]);
1921
1922           // Determine quadrant occupancies
1923           switch (op) {
1924             case GPC.Diff:
1925             case GPC.Int:
1926               tr = (inn[CLIP]) && (inn[SUBJ]);
1927               tl = (inn[CLIP] ^ e1.bundle[ABOVE][CLIP]) && (inn[SUBJ] ^ e1.bundle[ABOVE][SUBJ]);
1928               br = (inn[CLIP] ^ e0.bundle[ABOVE][CLIP]) && (inn[SUBJ] ^ e0.bundle[ABOVE][SUBJ]);
1929               bl = (inn[CLIP] ^ e1.bundle[ABOVE][CLIP] ^ e0.bundle[ABOVE][CLIP]) && (inn[SUBJ] ^ e1.bundle[ABOVE][SUBJ] ^ e0.bundle[ABOVE][SUBJ]);
1930               break;
1931             case GPC.Xor:
1932               tr = (inn[CLIP]) ^ (inn[SUBJ]);
1933               tl = (inn[CLIP] ^ e1.bundle[ABOVE][CLIP]) ^ (inn[SUBJ] ^ e1.bundle[ABOVE][SUBJ]);
1934               br = (inn[CLIP] ^ e0.bundle[ABOVE][CLIP]) ^ (inn[SUBJ] ^ e0.bundle[ABOVE][SUBJ]);
1935               bl = (inn[CLIP] ^ e1.bundle[ABOVE][CLIP] ^ e0.bundle[ABOVE][CLIP]) ^ (inn[SUBJ] ^ e1.bundle[ABOVE][SUBJ] ^ e0.bundle[ABOVE][SUBJ]);
1936               break;
1937             case GPC.Union:
1938               tr = (inn[CLIP]) || (inn[SUBJ]);
1939               tl = (inn[CLIP] ^ e1.bundle[ABOVE][CLIP]) || (inn[SUBJ] ^ e1.bundle[ABOVE][SUBJ]);
1940               br = (inn[CLIP] ^ e0.bundle[ABOVE][CLIP]) || (inn[SUBJ] ^ e0.bundle[ABOVE][SUBJ]);
1941               bl = (inn[CLIP] ^ e1.bundle[ABOVE][CLIP] ^ e0.bundle[ABOVE][CLIP]) || (inn[SUBJ] ^ e1.bundle[ABOVE][SUBJ] ^ e0.bundle[ABOVE][SUBJ]);
1942               break;
1943             default: assert(0);
1944           }
1945
1946           vclass = tr + (tl << 1) + (br << 2) + (bl << 3);
1947
1948           switch (vclass) {
1949             case EMN:
1950               new_tristrip(&tlist, e1, ix, iy);
1951               e0.outp[ABOVE] = e1.outp[ABOVE];
1952               break;
1953             case ERI:
1954               if (p) {
1955                 P_EDGE(prev_edge, e0, ABOVE, px, iy);
1956                 VERTEX(prev_edge, ABOVE, LEFT, px, iy);
1957                 VERTEX(e0, ABOVE, RIGHT, ix, iy);
1958                 e1.outp[ABOVE] = e0.outp[ABOVE];
1959                 e0.outp[ABOVE] = null;
1960               }
1961               break;
1962             case ELI:
1963               if (q) {
1964                 N_EDGE(next_edge, e1, ABOVE, nx, iy);
1965                 VERTEX(e1, ABOVE, LEFT, ix, iy);
1966                 VERTEX(next_edge, ABOVE, RIGHT, nx, iy);
1967                 e0.outp[ABOVE] = e1.outp[ABOVE];
1968                 e1.outp[ABOVE] = null;
1969               }
1970               break;
1971             case EMX:
1972               if (p && q) {
1973                 VERTEX(e0, ABOVE, LEFT, ix, iy);
1974                 e0.outp[ABOVE] = null;
1975                 e1.outp[ABOVE] = null;
1976               }
1977               break;
1978             case IMN:
1979               P_EDGE(prev_edge, e0, ABOVE, px, iy);
1980               VERTEX(prev_edge, ABOVE, LEFT, px, iy);
1981               N_EDGE(next_edge, e1, ABOVE, nx, iy);
1982               VERTEX(next_edge, ABOVE, RIGHT, nx, iy);
1983               new_tristrip(&tlist, prev_edge, px, iy);
1984               e1.outp[ABOVE] = prev_edge.outp[ABOVE];
1985               VERTEX(e1, ABOVE, RIGHT, ix, iy);
1986               new_tristrip(&tlist, e0, ix, iy);
1987               next_edge.outp[ABOVE] = e0.outp[ABOVE];
1988               VERTEX(next_edge, ABOVE, RIGHT, nx, iy);
1989               break;
1990             case ILI:
1991               if (p) {
1992                 VERTEX(e0, ABOVE, LEFT, ix, iy);
1993                 N_EDGE(next_edge, e1, ABOVE, nx, iy);
1994                 VERTEX(next_edge, ABOVE, RIGHT, nx, iy);
1995                 e1.outp[ABOVE] = e0.outp[ABOVE];
1996                 e0.outp[ABOVE] = null;
1997               }
1998               break;
1999             case IRI:
2000               if (q) {
2001                 VERTEX(e1, ABOVE, RIGHT, ix, iy);
2002                 P_EDGE(prev_edge, e0, ABOVE, px, iy);
2003                 VERTEX(prev_edge, ABOVE, LEFT, px, iy);
2004                 e0.outp[ABOVE] = e1.outp[ABOVE];
2005                 e1.outp[ABOVE] = null;
2006               }
2007               break;
2008             case IMX:
2009               if (p && q) {
2010                 VERTEX(e0, ABOVE, RIGHT, ix, iy);
2011                 VERTEX(e1, ABOVE, LEFT, ix, iy);
2012                 e0.outp[ABOVE] = null;
2013                 e1.outp[ABOVE] = null;
2014                 P_EDGE(prev_edge, e0, ABOVE, px, iy);
2015                 VERTEX(prev_edge, ABOVE, LEFT, px, iy);
2016                 new_tristrip(&tlist, prev_edge, px, iy);
2017                 N_EDGE(next_edge, e1, ABOVE, nx, iy);
2018                 VERTEX(next_edge, ABOVE, RIGHT, nx, iy);
2019                 next_edge.outp[ABOVE] = prev_edge.outp[ABOVE];
2020                 VERTEX(next_edge, ABOVE, RIGHT, nx, iy);
2021               }
2022               break;
2023             case IMM:
2024               if (p && q) {
2025                 VERTEX(e0, ABOVE, RIGHT, ix, iy);
2026                 VERTEX(e1, ABOVE, LEFT, ix, iy);
2027                 P_EDGE(prev_edge, e0, ABOVE, px, iy);
2028                 VERTEX(prev_edge, ABOVE, LEFT, px, iy);
2029                 new_tristrip(&tlist, prev_edge, px, iy);
2030                 N_EDGE(next_edge, e1, ABOVE, nx, iy);
2031                 VERTEX(next_edge, ABOVE, RIGHT, nx, iy);
2032                 e1.outp[ABOVE] = prev_edge.outp[ABOVE];
2033                 VERTEX(e1, ABOVE, RIGHT, ix, iy);
2034                 new_tristrip(&tlist, e0, ix, iy);
2035                 next_edge.outp[ABOVE] = e0.outp[ABOVE];
2036                 VERTEX(next_edge, ABOVE, RIGHT, nx, iy);
2037               }
2038               break;
2039             case EMM:
2040               if (p && q) {
2041                 VERTEX(e0, ABOVE, LEFT, ix, iy);
2042                 new_tristrip(&tlist, e1, ix, iy);
2043                 e0.outp[ABOVE]= e1.outp[ABOVE];
2044               }
2045               break;
2046             default:
2047               break;
2048           } // End of switch
2049         } // End of contributing intersection conditional
2050
2051         // Swap bundle sides in response to edge crossing
2052         if (e0.bundle[ABOVE][CLIP]) e1.bside[CLIP] = !e1.bside[CLIP];
2053         if (e1.bundle[ABOVE][CLIP]) e0.bside[CLIP] = !e0.bside[CLIP];
2054         if (e0.bundle[ABOVE][SUBJ]) e1.bside[SUBJ] = !e1.bside[SUBJ];
2055         if (e1.bundle[ABOVE][SUBJ]) e0.bside[SUBJ] = !e0.bside[SUBJ];
2056
2057         // Swap the edge bundles in the aet
2058         swap_intersecting_edge_bundles(&aet, intersect);
2059       } // End of IT loop
2060
2061       // Prepare for next scanbeam
2062       for (edge = aet; edge !is null; edge = next_edge) {
2063         next_edge = edge.next;
2064         succ_edge = edge.succ;
2065
2066         if (edge.top.y == yt && succ_edge) {
2067           // Replace AET edge by its successor
2068           succ_edge.outp[BELOW] = edge.outp[ABOVE];
2069           succ_edge.bstate[BELOW] = edge.bstate[ABOVE];
2070           succ_edge.bundle[BELOW][CLIP] = edge.bundle[ABOVE][CLIP];
2071           succ_edge.bundle[BELOW][SUBJ] = edge.bundle[ABOVE][SUBJ];
2072           prev_edge = edge.prev;
2073           if (prev_edge) prev_edge.next = succ_edge; else aet = succ_edge;
2074           if (next_edge) next_edge.prev = succ_edge;
2075           succ_edge.prev = prev_edge;
2076           succ_edge.next = next_edge;
2077         } else {
2078           // Update this edge
2079           edge.outp[BELOW] = edge.outp[ABOVE];
2080           edge.bstate[BELOW] = edge.bstate[ABOVE];
2081           edge.bundle[BELOW][CLIP] = edge.bundle[ABOVE][CLIP];
2082           edge.bundle[BELOW][SUBJ] = edge.bundle[ABOVE][SUBJ];
2083           edge.xb = edge.xt;
2084         }
2085         edge.outp[ABOVE] = null;
2086       }
2087     }
2088   } // === END OF SCANBEAM PROCESSING ==================================
2089
2090   // Generate result tristrip from tlist
2091   result.strip = null;
2092   result.num_strips = count_tristrips(tlist);
2093   if (result.num_strips > 0) {
2094     result.strip = MALLOC!gpc_vertex_list(result.num_strips);
2095     s = 0;
2096     for (tn = tlist; tn !is null; tn = tnn) {
2097       tnn = tn.next;
2098       if (tn.active > 2) {
2099         // Valid tristrip: copy the vertices and free the heap
2100         result.strip[s].num_vertices= tn.active;
2101         result.strip[s].vertex = MALLOC!gpc_vertex(tn.active);
2102         v = 0;
2103         if (GPC_INVERT_TRISTRIPS) {
2104           lt = tn.v[RIGHT];
2105           rt = tn.v[LEFT];
2106         } else {
2107           lt = tn.v[LEFT];
2108           rt = tn.v[RIGHT];
2109         }
2110         while (lt || rt) {
2111           if (lt) {
2112             ltn = lt.next;
2113             result.strip[s].vertex[v].x = lt.x;
2114             result.strip[s].vertex[v].y = lt.y;
2115             ++v;
2116             FREE(lt);
2117             lt = ltn;
2118           }
2119           if (rt) {
2120             rtn = rt.next;
2121             result.strip[s].vertex[v].x = rt.x;
2122             result.strip[s].vertex[v].y = rt.y;
2123             ++v;
2124             FREE(rt);
2125             rt = rtn;
2126           }
2127         }
2128         ++s;
2129       } else {
2130         // Invalid tristrip: just free the heap
2131         for (lt = tn.v[LEFT]; lt !is null; lt = ltn) {
2132           ltn = lt.next;
2133           FREE(lt);
2134         }
2135         for (rt = tn.v[RIGHT]; rt !is null; rt = rtn) {
2136           rtn = rt.next;
2137           FREE(rt);
2138         }
2139       }
2140       FREE(tn);
2141     }
2142   }
2143
2144   // Tidy up
2145   reset_it(&it);
2146   reset_lmt(&lmt);
2147   FREE(c_heap);
2148   FREE(s_heap);
2149   FREE(sbt);
2150 }