/*******************************************************************************
*                                                                              *
* Author: Alexis Naveros                                                       *
* Algorithm and implementation inspired by the Clipper2 library                *
* This is partly a rewrite of Clipper2, more optimized, scalable and robust    *
*                                                                              *
* Clipper2's Copyright: (C) Angus Johnson 2010-2022                            *
* Clipper2's copyright notice included below                                   *
*                                                                              *
*******************************************************************************/
/*******************************************************************************
* Author    :  Angus Johnson                                                   *
* Version   :  Clipper2 - beta                                                 *
* Website   :  http://www.angusj.com                                           *
* Copyright :  Angus Johnson 2010-2022                                         *
* License   :  http://www.boost.org/LICENSE_1_0.txt                            *
*******************************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <stddef.h>
#include <stdint.h>
#include <string.h>
#include <math.h>
#include <float.h>

#include "cpusimd.h"

#include "mmslab.h"
#include "mmheap.h"

#include "mathlineline.h"
#include "mathshewchuk.h"

#include "polyclip.h"
#include "polywithin.h"


////


/* Improve cache locality, reduce overhead, a good idea all around */
#define CLP_CONFIG_ENABLE_SLAB_ALLOCATION (1)

/* A good idea when the count of items to sort gets high */
#define CLP_CONFIG_ENABLE_RADIX_SORT (1)

/* What's the count threshold to switch from hybrid sort to radix sort? */
#define CLP_CONFIG_RADIX_SORT_COUNT_THRESHOLD (4096)

/* Do we want to also sort intersection nodes on X? It shouldn't alter the validity of the results */
#define CLP_CONFIG_SORT_INODE_X (0)

/* We allocate extra vertices beyond the "vertexcount" range in order to allow guard values rather than bound checks */
#define CLP_CONFIG_PATHVLIST_EXTRA_ALLOC (PW_STORAGE_PADDING_BYTES/(2*8))


#if defined(__GNUC__) || defined(__INTEL_COMPILER) || defined(__clang__)
 #define CC_ALWAYSINLINE __attribute__((always_inline))
#else
 #define CC_ALWAYSINLINE
#endif


////


#define CLP_EPSILON (1.0e-12)


typedef struct
{
  int64_t x;
  int64_t y;
} clpPoint64;

typedef struct
{
  double x;
  double y;
} clpPointD;

static inline clpPoint64 clpBuildPoint64( int64_t x, int64_t y )
{
  clpPoint64 p64;
  p64.x = x;
  p64.y = y;
  return p64;
}

static inline clpPointD clpBuildPointD( double x, double y )
{
  clpPointD pd;
  pd.x = x;
  pd.y = y;
  return pd;
}


typedef struct clpActiveEdge clpActiveEdge;
typedef struct clpVertex clpVertex;
typedef struct clpLocalMinima clpLocalMinima;
typedef struct clpOutRec clpOutRec;
typedef struct clpJoiner clpJoiner;

typedef struct
{
  clpPoint64 pt;
  clpActiveEdge *edge1;
  clpActiveEdge *edge2;
} clpIntersectNode;



#define MM_VECTOR_ITEMTYPE clpLocalMinima *
#define MM_VECTOR_TYPENAME clpLocalMinimaVector
#define MM_VECTOR_PREFIX clpLocalMinimaVector_
#define MM_VECTOR_OPTION_ORDERED (1)
#define MM_VECTOR_OPTION_DEBUG_CALLS (0)
#include "mmvector.h"

#define MM_VECTOR_ITEMTYPE clpOutRec *
#define MM_VECTOR_TYPENAME clpOutRecVector
#define MM_VECTOR_PREFIX clpOutRecVector_
#define MM_VECTOR_OPTION_ORDERED (1)
#define MM_VECTOR_OPTION_DEBUG_CALLS (0)
#include "mmvector.h"

#define MM_VECTOR_ITEMTYPE clpIntersectNode
#define MM_VECTOR_TYPENAME clpIntersectNodeVector
#define MM_VECTOR_PREFIX clpIntersectNodeVector_
#define MM_VECTOR_OPTION_ORDERED (1)
#define MM_VECTOR_OPTION_DEBUG_CALLS (0)
#include "mmvector.h"

#define MM_VECTOR_ITEMTYPE clpJoiner *
#define MM_VECTOR_TYPENAME clpJoinerVector
#define MM_VECTOR_PREFIX clpJoinerVector_
#define MM_VECTOR_OPTION_ORDERED (1)
#define MM_VECTOR_OPTION_DEBUG_CALLS (0)
#include "mmvector.h"

#define MM_VECTOR_ITEMTYPE clpVertex *
#define MM_VECTOR_TYPENAME clpVertexList
#define MM_VECTOR_PREFIX clpVertexList_
#define MM_VECTOR_OPTION_ORDERED (1)
#define MM_VECTOR_OPTION_DEBUG_CALLS (0)
#include "mmvector.h"


#define CLP_VERTEX_FLAG_OPENSTART (0x1)
#define CLP_VERTEX_FLAG_OPENEND (0x2)
#define CLP_VERTEX_FLAG_LOCALMAX (0x4)
#define CLP_VERTEX_FLAG_LOCALMIN (0x8)

struct clpVertex
{
  clpPoint64 pt;
  clpVertex *next;
  clpVertex *prev;
  uint32_t flags;
  uint32_t unused; /* Wasteful implicit padding... let's pack bits in the first bits of the pointers? */
};

typedef struct clpOutPt clpOutPt;

struct clpOutPt
{
  clpPoint64 pt;
  clpOutPt *next;
  clpOutPt *prev;
  clpOutRec *outrec;
  clpJoiner *joiner;
};


/* clpOutRec: contains a path in the clipping solution. Edges in the AEL will */
/* have clpOutRec pointers assigned when they form part of the clipping solution */
struct clpOutRec
{
  size_t idx;
  clpOutRec *owner;
  clpOutRecVector splits;
  clpActiveEdge *front_edge;
  clpActiveEdge *back_edge;
  clpOutPt *pts;
  clpTree *shapetree;
  int64_t bbox[4];
  clpPath path;
  uint8_t openflag;
  /* 7 bytes of wasted alignment storage here */
};


struct clpActiveEdge
{
  clpPoint64 bot;
  clpPoint64 top;
  int64_t scanbeamx;    /* Current (updated at every new scanline) ~ TODO: maintain inactive queue */
  double slopedx;
  uint8_t leftboundflag;
  int8_t windingdelta;      /* 1 or -1 depending on winding direction */
  /* 6 bytes of wasteful padding here */
  int windingcount;
  int windingcount2;    /* winding count of the opposite polytype */
  clpOutRec *outrec;

  /* Linked list of all edges (from left to right) that are present */
  /* (or 'active') within the current scanbeam (a horizontal 'beam' that */
  /* sweeps from bottom to top over the paths in the clipping operation). */
  clpActiveEdge *aelprev;
  clpActiveEdge *aelnext;

  /* Linked list used when sorting edges into their new positions at the */
  /* top of scanbeams, but also (re)used to process horizontals. */
  clpActiveEdge *selprev;
  clpActiveEdge *selnext;

  clpActiveEdge *jump;
  clpVertex *vertextop;
  clpLocalMinima *local_min;  /* the bottom of an edge 'bound' */
};

struct clpLocalMinima
{
  clpVertex *vertex;
  uint8_t polytype;
  uint8_t openflag;
};

struct clpJoiner
{
  int idx;
  clpOutPt *op1;
  clpOutPt *op2;
  clpJoiner *next1;
  clpJoiner *next2;
  clpJoiner *nextH;
};


////


struct clpClipper
{
  /* Runtime status */
  int64_t scanbeambottomy;
  uint8_t cliptype;
  uint8_t fillrule;
  uint8_t fillpos;
  uint8_t openpathsflag;
  uint8_t minimalistsorted;
  uint8_t polytreeflag;
  uint8_t succeessflag;

#if CLP_CONFIG_ENABLE_SLAB_ALLOCATION
  mmSlab slaboutpt;
  mmSlab slablm;
  mmSlab slaboutrec;
  mmSlab slabjoiner;
  mmSlab slabactive;
#endif

  /* Linked list of all edges (from left to right) that are present */
  /* (or 'active') within the current scanbeam (a horizontal 'beam' that */
  /* sweeps from bottom to top over the paths in the clipping operation). */
  clpActiveEdge *activeedgelist;

  /* Linked list used when sorting edges into their new positions at the */
  /* top of scanbeams, but also (re)used to process horizontals. */
  clpActiveEdge *sortededgelist;
  clpJoiner *horzjoinerlist;

  /* Runtime data storage */
  clpVertexList vertexlist;
  clpLocalMinimaVector minimalist;
  size_t minimalistwalkindex;
  mmHeap64sb scanlinelist;
  clpIntersectNodeVector intersectnodelist;
  clpOutRecVector outreclist;
  clpJoinerVector joinerlist;

  /* Options to expose through the API */
  uint8_t orientationreversedflag;
  uint8_t preservecollinearflag;
  uint8_t reversesolutionflag;
};


////


static inline void DebugDebugPrintPtChain( clpOutPt *op, const char *header )
{
  clpOutPt *startop;
  startop = op;
  printf( "%s <<<<<<\n", header );
  for( ; ; )
  {
    printf( "%s %d %d\n", header, (int)op->pt.x, (int)op->pt.y );
    if( op->joiner )
      break;
    op = op->next;
    if( op == startop )
      break;
  }
  return;
}

static inline void DebugDebugPrintIntersectList( clpClipper *clp, const char *header )
{
  clpActiveEdge *edge;
  printf( "%s <<<<<<\n", header );
  size_t nodeindex, nodecount;
  clpIntersectNode *node0;
  nodecount = clpIntersectNodeVector_GetCount( &clp->intersectnodelist );
  for( nodeindex = 0 ; nodeindex < nodecount ; nodeindex++ )
  {
    node0 = clpIntersectNodeVector_Access( &clp->intersectnodelist, nodeindex );
    printf( "%s %d %d\n", header, (int)node0->pt.x, (int)node0->pt.y );
    edge = node0->edge1;
    printf( "  bot %d %d ~ top %d %d ~ curr_x %d ~ dx %f\n", (int)edge->bot.x, (int)edge->bot.y, (int)edge->top.x, (int)edge->top.y, (int)edge->scanbeamx, edge->slopedx );
    edge = node0->edge2;
    printf( "  bot %d %d ~ top %d %d ~ curr_x %d ~ dx %f\n", (int)edge->bot.x, (int)edge->bot.y, (int)edge->top.x, (int)edge->top.y, (int)edge->scanbeamx, edge->slopedx );
  }
  return;
}


////


static inline int clpPoint64Equal( const clpPoint64 *a, const clpPoint64 *b )
{
  return ( ( ( a->x - b->x ) | ( a->y - b->y ) ) == 0 );
}

static inline double clpCrossProduct( const clpPoint64 *pt1, const clpPoint64 *pt2, const clpPoint64 *pt3 )
{
  return ( ( (double)( pt2->x - pt1->x ) * (double)( pt3->y - pt2->y ) ) - ( (double)( pt2->y - pt1->y ) * (double)( pt3->x - pt2->x ) ) );
}

static inline double clpDotProduct( const clpPoint64 *pt1, const clpPoint64 *pt2, const clpPoint64 *pt3 )
{
  return ( ( (double)( pt2->x - pt1->x ) * (double)( pt3->x - pt2->x ) ) + ( (double)( pt2->y - pt1->y ) * (double)( pt3->y - pt2->y ) ) );
}

static inline double clpDistanceSqr( const clpPoint64 *pt1, const clpPoint64 *pt2 )
{
  double dx, dy;
  dx = (double)( pt1->x - pt2->x );
  dy = (double)( pt1->y - pt2->y );
#if CLP_DEBUG_VERBOSE
  printf( "DistanceSqr %f\n", ( dx * dx ) + ( dy * dy ) );
#endif
  return ( dx * dx ) + ( dy * dy );
}

static inline void clpBboxSetEmpty( int64_t *bbox )
{
  int64_t bbminx, bbmaxx, bbminy, bbmaxy;
  bbminx = INT64_MAX;
  bbmaxx = INT64_MIN;
  bbminy = INT64_MAX;
  bbmaxy = INT64_MIN;
  bbox[0] = bbminx;
  bbox[1] = bbmaxx;
  bbox[2] = bbminy;
  bbox[3] = bbmaxy;
  return;
}


////


static inline clpOutPt *clpNewOutPt( clpClipper *clp, const clpPoint64 *pt, clpOutRec *outrec )
{
  clpOutPt *outpt;
#if CLP_CONFIG_ENABLE_SLAB_ALLOCATION
  outpt = mmSlabAlloc( &clp->slaboutpt );
#else
  outpt = (clpOutPt *)malloc( sizeof(clpOutPt) );
#endif
  outpt->pt = *pt;
  outpt->next = outpt;
  outpt->prev = outpt;
  outpt->outrec = outrec;
  outpt->joiner = 0;
  return outpt;
}

static inline void clpFreeOutPt( clpClipper *clp, clpOutPt *outpt )
{
#if CLP_CONFIG_ENABLE_SLAB_ALLOCATION
  mmSlabRelease( &clp->slaboutpt, outpt );
#else
  free( (void *)outpt );
#endif
  return;
}


static inline clpLocalMinima *clpNewLocalMinima( clpClipper *clp, clpVertex *vertex, uint8_t polytype, uint8_t openflag )
{
  clpLocalMinima *lm;
#if CLP_CONFIG_ENABLE_SLAB_ALLOCATION
  lm = mmSlabAlloc( &clp->slablm );
#else
  lm = (clpLocalMinima *)malloc( sizeof(clpLocalMinima) );
#endif
  lm->vertex = vertex;
  lm->polytype = polytype;
  lm->openflag = openflag;
  return lm;
}

static inline void clpFreeLocalMinima( clpClipper *clp, clpLocalMinima *lm )
{
#if CLP_CONFIG_ENABLE_SLAB_ALLOCATION
  mmSlabRelease( &clp->slablm, lm );
#else
  free( (void *)lm );
#endif
  return;
}


static inline clpOutRec *clpNewOutRec( clpClipper *clp )
{
  clpOutRec *outrec;
#if CLP_CONFIG_ENABLE_SLAB_ALLOCATION
  outrec = mmSlabAlloc( &clp->slaboutrec );
#else
  outrec = (clpOutRec *)malloc( sizeof(clpOutRec) );
#endif
  memset( (void *)outrec, 0, sizeof(clpOutRec) );
  clpOutRecVector_Init( &outrec->splits, 0 );
  clpBboxSetEmpty( outrec->bbox );
  return outrec;
}

static inline void clpFreeOutRec( clpClipper *clp, clpOutRec *outrec )
{
  clpOutRecVector_Clear( &outrec->splits );
#if CLP_CONFIG_ENABLE_SLAB_ALLOCATION
  mmSlabRelease( &clp->slaboutrec, outrec );
#else
  free( (void *)outrec );
#endif
  return;
}


static inline clpJoiner *clpNewJoiner( clpClipper *clp, clpOutPt *op1, clpOutPt *op2, clpJoiner* nexth )
{
  clpJoiner *joiner;
#if CLP_CONFIG_ENABLE_SLAB_ALLOCATION
  joiner = mmSlabAlloc( &clp->slabjoiner );
#else
  joiner = (clpJoiner *)malloc( sizeof(clpJoiner) );
#endif
  joiner->idx = -1;
  joiner->op1 = op1;
  joiner->op2 = op2;
  joiner->nextH = nexth;
  joiner->next1 = op1->joiner;
  op1->joiner = joiner;
  if( op2 )
  {
    joiner->next2 = op2->joiner;
    op2->joiner = joiner;
  }
  else
    joiner->next2 = 0;
  return joiner;
}

static inline void clpFreeJoiner( clpClipper *clp, clpJoiner *joiner )
{
#if CLP_CONFIG_ENABLE_SLAB_ALLOCATION
  mmSlabRelease( &clp->slabjoiner, joiner );
#else
  free( (void *)joiner );
#endif
  return;
}


static inline clpActiveEdge *clpNewActiveEdge( clpClipper *clp )
{
  clpActiveEdge *active;
  /* TODO: Replace by array allocation, fixed indices, unordered array ? */
  /* --> Replace all storage pointers to clpActiveEdge with uint32_t indices into fixed unordered array <-- */
#if CLP_CONFIG_ENABLE_SLAB_ALLOCATION
  active = mmSlabAlloc( &clp->slabactive );
#else
  active = (clpActiveEdge *)malloc( sizeof(clpActiveEdge) );
#endif
  memset( (void *)active, 0, sizeof(clpActiveEdge) );
  active->scanbeamx = 0;    /* current (updated at every new scanline) ~ TODO: inactive queue! */
  active->slopedx = 0.0;
  active->windingdelta = 1;     /* 1 or -1 depending on winding direction */
  active->windingcount = 0;
  active->windingcount2 = 0;    /* winding count of the opposite polytype */
  active->leftboundflag = 0;
  return active;
}

static inline void clpFreeActiveEdge( clpClipper *clp, clpActiveEdge *active )
{
#if CLP_CONFIG_ENABLE_SLAB_ALLOCATION
  mmSlabRelease( &clp->slabactive, active );
#else
  free( (void *)active );
#endif
  return;
}


////


/* Every closed path (or polygon) is made up of a series of vertices forming */
/* edges that alternate between going up (relative to the Y-axis) and going */
/* down. Edges consecutively going up or consecutively going down are called */
/* 'bounds' (or sides if they're simple polygons). 'Local Minima' refer to */
/* vertices where descending bounds become ascending ones. */


#define RSORT_MAIN clpLocalMinimaRadixSortX
#define RSORT_RADIX clpLocalMinimaRadixSortXRadix
#define RSORT_TYPE clpLocalMinima *
#define RSORT_RADIXBITS (11)
#define RSORT_BIGGEST_FIRST (1)
#define RSORT_TESTFULLBIN (1)
static inline CC_ALWAYSINLINE uint32_t clpLocalMinimaRadixSortXRadix( clpLocalMinima **lmp, int bitindex )
{
  uint32_t radix;
  uint64_t sortx;
  clpLocalMinima *lm;
  lm = *lmp;
  /* TODO: We should do that XOR before and after the sort, not in here */
  sortx = (uint64_t)(lm->vertex)->pt.x ^ 0x8000000000000000;
  radix = ( sortx >> bitindex );
  return radix & ((1<<RSORT_RADIXBITS)-1);
}
#include "ccradixsort.h"
#undef RSORT_MAIN
#undef RSORT_RADIX
#undef RSORT_TYPE
#undef RSORT_RADIXBITS
#undef RSORT_BIGGEST_FIRST
#undef RSORT_TESTFULLBIN


#define RSORT_MAIN clpLocalMinimaRadixSortY
#define RSORT_RADIX clpLocalMinimaRadixSortYRadix
#define RSORT_TYPE clpLocalMinima *
#define RSORT_RADIXBITS (11)
#define RSORT_BIGGEST_FIRST (1)
#define RSORT_TESTFULLBIN (1)
static inline CC_ALWAYSINLINE uint32_t clpLocalMinimaRadixSortYRadix( clpLocalMinima **lmp, int bitindex )
{
  uint32_t radix;
  uint64_t sorty;
  clpLocalMinima *lm;
  lm = *lmp;
  /* TODO: We should do that XOR before and after the sort, not in here */
  sorty = (uint64_t)(lm->vertex)->pt.y ^ 0x8000000000000000;
  radix = ( sorty >> bitindex );
  return radix & ((1<<RSORT_RADIXBITS)-1);
}
#include "ccradixsort.h"
#undef RSORT_MAIN
#undef RSORT_RADIX
#undef RSORT_TYPE
#undef RSORT_RADIXBITS
#undef RSORT_BIGGEST_FIRST
#undef RSORT_TESTFULLBIN


#define HSORT_MAIN clpLocalMinimaHybridSort
#define HSORT_CMP clpLocalMinimaHybridSortCmp
#define HSORT_TYPE clpLocalMinima *
static inline CC_ALWAYSINLINE int clpLocalMinimaHybridSortCmp( clpLocalMinima **v0, clpLocalMinima **v1 )
{
  clpVertex *vertex0, *vertex1;
  vertex0 = (*v0)->vertex;
  vertex1 = (*v1)->vertex;
  return ( ( vertex0->pt.y != vertex1->pt.y ) ? ( vertex0->pt.y < vertex1->pt.y ) : ( vertex0->pt.x < vertex1->pt.x ) );
}
#if CPU_SSE4_2_SUPPORT
 #define HSORT_CMP_XMM clpLocalMinimaHybridSortCmpXmm
static inline CC_ALWAYSINLINE __m128 clpLocalMinimaHybridSortCmpXmm( clpLocalMinima **v0, clpLocalMinima **v1 )
{
  clpVertex *vertex0, *vertex1;
  __m128i v0x, v0y, v1x, v1y;
  __m128i vcmp;
  vertex0 = (*v0)->vertex;
  vertex1 = (*v1)->vertex;
  v0x = _mm_loadl_epi64( (void *)&vertex0->pt.x );
  v0y = _mm_loadl_epi64( (void *)&vertex0->pt.y );
  v1x = _mm_loadl_epi64( (void *)&vertex1->pt.x );
  v1y = _mm_loadl_epi64( (void *)&vertex1->pt.y );
  vcmp = _mm_cmpgt_epi64( v0y, v1y );
  vcmp = _mm_or_si128( _mm_cmpeq_epi64( v0y, v1y ), _mm_cmpgt_epi64( v0x, v1x ) );
  return _mm_castsi128_ps( vcmp );
}
#endif
#include "cchybridsort.h"
#undef HSORT_MAIN
#undef HSORT_CMP
#undef HSORT_CMP_XMM
#undef HSORT_TYPE


#define RSORT_MAIN clpIntersectNodeRadixSortX
#define RSORT_RADIX clpIntersectNodeRadixSortXRadix
#define RSORT_TYPE clpIntersectNode
#define RSORT_RADIXBITS (11)
#define RSORT_BIGGEST_FIRST (1)
#define RSORT_TESTFULLBIN (1)
static inline CC_ALWAYSINLINE uint32_t clpIntersectNodeRadixSortXRadix( clpIntersectNode *node, int bitindex )
{
  uint32_t radix;
  uint64_t sortx;
  /* TODO: We should do that XOR before and after the sort, not in here */
  sortx = (uint64_t)node->pt.x ^ 0x8000000000000000;
  radix = ( sortx >> bitindex );
  return radix & ((1<<RSORT_RADIXBITS)-1);
}
#include "ccradixsort.h"
#undef RSORT_MAIN
#undef RSORT_RADIX
#undef RSORT_TYPE
#undef RSORT_RADIXBITS
#undef RSORT_BIGGEST_FIRST
#undef RSORT_TESTFULLBIN


#define RSORT_MAIN clpIntersectNodeRadixSortY
#define RSORT_RADIX clpIntersectNodeRadixSortYRadix
#define RSORT_TYPE clpIntersectNode
#define RSORT_RADIXBITS (11)
#define RSORT_BIGGEST_FIRST (1)
#define RSORT_TESTFULLBIN (1)
static inline CC_ALWAYSINLINE uint32_t clpIntersectNodeRadixSortYRadix( clpIntersectNode *node, int bitindex )
{
  uint32_t radix;
  uint64_t sorty;
  /* TODO: We should do that XOR before and after the sort, not in here */
  sorty = (uint64_t)node->pt.y ^ 0x8000000000000000;
  radix = ( sorty >> bitindex );
  return radix & ((1<<RSORT_RADIXBITS)-1);
}
#include "ccradixsort.h"
#undef RSORT_MAIN
#undef RSORT_RADIX
#undef RSORT_TYPE
#undef RSORT_RADIXBITS
#undef RSORT_BIGGEST_FIRST
#undef RSORT_TESTFULLBIN


#define HSORT_MAIN clpIntersectNodeHybridSort
#define HSORT_CMP clpIntersectNodeHybridSortCmp
#define HSORT_TYPE clpIntersectNode
#define HSORT_UNROLL (1)
#define HSORT_USE_TSORT (1)
static inline CC_ALWAYSINLINE int clpIntersectNodeHybridSortCmp( clpIntersectNode *node0, clpIntersectNode *node1 )
{
 #if CLP_CONFIG_SORT_INODE_X
  return ( ( node0->pt.y != node1->pt.y ) ? ( node0->pt.y < node1->pt.y ) : ( node0->pt.x > node1->pt.x ) );
 #else
  return ( node0->pt.y < node1->pt.y );
 #endif
}
#if CPU_SSE4_2_SUPPORT && !CLP_CONFIG_SORT_INODE_X
 #define HSORT_CMP_XMM clpIntersectNodeHybridSortCmpXmm
static inline CC_ALWAYSINLINE __m128 clpIntersectNodeHybridSortCmpXmm( clpIntersectNode *node0, clpIntersectNode *node1 )
{
  __m128i vcmp;
  vcmp = _mm_cmpgt_epi64( _mm_loadl_epi64( (void *)&node0->pt.y ), _mm_loadl_epi64( (void *)&node1->pt.y ) );
  return _mm_castsi128_ps( vcmp );
}
#endif
#include "cchybridsort.h"
#undef HSORT_MAIN
#undef HSORT_CMP
#undef HSORT_CMP_XMM
#undef HSORT_TYPE
#undef HSORT_UNROLL
#undef HSORT_USE_TSORT


////


static inline int clpIsHotEdge( const clpActiveEdge *e )
{
  return ( e->outrec != 0 );
}

static inline int clpIsActiveEdgeOpen( const clpActiveEdge *e )
{
  return ( (e->local_min)->openflag != 0 );
}

static inline uint32_t clpIsVertexOpenEnd( const clpVertex *v )
{
  return ( v->flags & ( CLP_VERTEX_FLAG_OPENSTART | CLP_VERTEX_FLAG_OPENEND ) );
}

static inline int clpIsActiveEdgeOpenEnd( const clpActiveEdge *ae )
{
  return clpIsVertexOpenEnd( ae->vertextop );
}

static inline clpActiveEdge *clpGetPrevHotEdge( const clpActiveEdge *e )
{
  clpActiveEdge *prev;
  prev = e->aelprev;
  while( prev && ( clpIsActiveEdgeOpen( prev ) || !clpIsHotEdge( prev ) ) )
    prev = prev->aelprev;
  return prev;
}

static inline int clpIsFront( const clpActiveEdge *e )
{
  return ( e == (e->outrec)->front_edge );
}


static inline int64_t clpGetActiveTopX( const clpActiveEdge *ae, const int64_t currentY )
{
  if( ( currentY == ae->top.y ) || ( ae->top.x == ae->bot.x ) )
    return ae->top.x;
  else if( currentY == ae->bot.y )
    return ae->bot.x;
  else
    return ae->bot.x + (int64_t)( round( ae->slopedx * (double)( currentY - ae->bot.y ) ) );
}


static inline int clpIsHorizontal( const clpActiveEdge *e )
{
  return ( e->top.y == e->bot.y );
}


static inline int clpIsHeadingRightHorz( const clpActiveEdge *e )
{
  return ( e->slopedx == -DBL_MAX );
}


static inline int clpIsHeadingLeftHorz( const clpActiveEdge *e )
{
  return ( e->slopedx == DBL_MAX );
}


static inline void clpSwapActives( clpActiveEdge **e1, clpActiveEdge **e2 )
{
  clpActiveEdge *e;
  e = *e1;
  *e1 = *e2;
  *e2 = e;
  return;
}


static inline uint8_t clpGetPolyType( const clpActiveEdge *e )
{
  return (e->local_min)->polytype;
}


static inline int clpIsSamePolyType( const clpActiveEdge *e1, const clpActiveEdge *e2 )
{
  return ( (e1->local_min)->polytype == (e2->local_min)->polytype );
}


static inline int clpGetIntersectPoint( const clpPoint64 *ln1a, const clpPoint64 *ln1b, const clpPoint64 *ln2a, const clpPoint64 *ln2b, clpPointD *ip )
{
  double m1, b1, m2, b2;

  if( ln1b->x == ln1a->x )
  {
    if( ln2b->x == ln2a->x )
    {
      ip->x = 0.0;
      ip->y = 0.0;
      return 0;
    }
    m2 = (double)( ln2b->y - ln2a->y ) / (double)( ln2b->x - ln2a->x );
    b2 = ln2a->y - m2 * ln2a->x;
    ip->x = (double)( ln1a->x );
    ip->y = m2 * ln1a->x + b2;
  }
  else if( ln2b->x == ln2a->x )
  {
    m1 = (double)( ln1b->y - ln1a->y ) / (double)( ln1b->x - ln1a->x );
    b1 = ln1a->y - m1 * ln1a->x;
    ip->x = (double)( ln2a->x );
    ip->y = m1 * ln2a->x + b1;
  }
  else
  {
    m1 = (double)( ln1b->y - ln1a->y ) / (double)( ln1b->x - ln1a->x );
    b1 = ln1a->y - m1 * ln1a->x;
    m2 = (double)( ln2b->y - ln2a->y ) / (double)( ln2b->x - ln2a->x );
    b2 = ln2a->y - m2 * ln2a->x;
    if( fabs( m1 - m2 ) > CLP_EPSILON )
    {
      ip->x = ( b2 - b1 ) / ( m1 - m2 );
      ip->y = m1 * ip->x + b1;
    }
    else
    {
      ip->x = (double)( ln1a->x + ln1b->x ) * 0.5;
      ip->y = (double)( ln1a->y + ln1b->y ) * 0.5;
    }
  }
  return 1;
}


static inline double clpGetSlopeDx( const clpPoint64 *pt1, const clpPoint64 *pt2 )
{
  int64_t dy;
  dy = pt2->y - pt1->y;
  if( dy != 0 )
    return (double)( pt2->x - pt1->x ) / (double)dy;
  else if( pt2->x > pt1->x )
    return -DBL_MAX;
  else
    return DBL_MAX;
}

static inline void clpActiveUpdateDx( clpActiveEdge *e )
{
  e->slopedx = clpGetSlopeDx( &e->bot, &e->top );
  return;
}


static inline clpVertex *clpGetNextVertex( const clpActiveEdge *ae )
{
  if( ae->windingdelta > 0 )
    return (ae->vertextop)->next;
  else
    return (ae->vertextop)->prev;
}


/* PrevPrevVertex: useful to get the (inverted Y-axis) top of the */
/* alternate edge (ie left or right bound) during edge insertion. */
static inline clpVertex *clpGetPrevPrevVertex( const clpActiveEdge *ae )
{
  if( ae->windingdelta > 0 )
    return ((ae->vertextop)->prev)->prev;
  else
    return ((ae->vertextop)->next)->next;
}


static inline clpActiveEdge *clpExtractFromSEL( clpActiveEdge *ae )
{
  clpActiveEdge *res;
  res = ae->selnext;
  if( res )
    res->selprev = ae->selprev;
  (ae->selprev)->selnext = res;
  return res;
}


static inline void clpInsert1Before2InSEL( clpActiveEdge *ae1, clpActiveEdge *ae2 )
{
  ae1->selprev = ae2->selprev;
  if( ae1->selprev )
    ae1->selprev->selnext = ae1;
  ae1->selnext = ae2;
  ae2->selprev = ae1;
  return;
}

static inline uint32_t clpIsVertexMaxima( const clpVertex *v )
{
  return ( v->flags & CLP_VERTEX_FLAG_LOCALMAX );
}


static inline int clpIsActiveEdgeMaxima( const clpActiveEdge *ae )
{
  return clpIsVertexMaxima( ae->vertextop );
}

static clpVertex *clpGetCurrYMaximaVertex( const clpActiveEdge *e )
{
  clpVertex *result;
  result = e->vertextop;
  if( e->windingdelta > 0 )
  {
    while( result->next->pt.y == result->pt.y )
      result = result->next;
  }
  else
  {
    while( result->prev->pt.y == result->pt.y )
      result = result->prev;
  }
  if( !clpIsVertexMaxima( result ) )
    result = 0; /* Not a maxima */
  return result;
}

static inline clpActiveEdge *clpGetMaximaPair( const clpActiveEdge *e )
{
  clpActiveEdge *e2;
  e2 = e->aelnext;
  while( e2 )
  {
    if( e2->vertextop == e->vertextop )
      return e2;
    e2 = e2->aelnext;
  }
  return 0;
}

static clpActiveEdge *clpGetHorzMaximaPair( const clpActiveEdge *horz, const clpVertex *vert_max )
{
  clpActiveEdge *result;
  /* We can't be sure whether the MaximaPair is on the left or right, so let's find out */
  result = horz->aelprev;
  while( result && ( result->scanbeamx >= vert_max->pt.x ) )
  {
    if( result->vertextop == vert_max )
      return result;
    result = result->aelprev;
  }
  result = horz->aelnext;
  while( result && ( clpGetActiveTopX( result, horz->top.y ) <= vert_max->pt.x ) )
  {
    if( result->vertextop == vert_max )
      return result;
    result = result->aelnext;
  }
  return 0;
}

static clpOutPt *clpInsertOp( clpClipper *clp, const clpPoint64 *pt, clpOutPt *insertafter )
{
  clpOutPt *result;
  result = clpNewOutPt( clp, pt, insertafter->outrec );
  result->next = insertafter->next;
  (insertafter->next)->prev = result;
  insertafter->next = result;
  result->prev = insertafter;
  return result;
}


static inline clpOutPt *clpDisposeOutPt( clpClipper *clp, clpOutPt *op )
{
  clpOutPt *result;
  result = op->next;
  (op->prev)->next = op->next;
  (op->next)->prev = op->prev;
  clpFreeOutPt( clp, op );
  return result;
}


static inline void clpDisposeOutPts( clpClipper *clp, clpOutRec *outrec )
{
  clpOutPt *op2, *tmp;
  if( !outrec->pts )
    return;
  op2 = outrec->pts->next;
  while( op2 != outrec->pts )
  {
    tmp = op2->next;
    clpFreeOutPt( clp, op2 );
    op2 = tmp;
  }
  clpFreeOutPt( clp, outrec->pts );
  outrec->pts = 0;
  return;
}


static inline void clpSetSides( clpOutRec *outrec, clpActiveEdge *start_edge, clpActiveEdge *end_edge )
{
  outrec->front_edge = start_edge;
  outrec->back_edge = end_edge;
  return;
}


static void clpSwapOutrecs( clpActiveEdge *e1, clpActiveEdge *e2 )
{
  clpActiveEdge *e;
  clpOutRec *or1, *or2;

  or1 = e1->outrec;
  or2 = e2->outrec;
  if( or1 == or2 )
  {
    e = or1->front_edge;
    or1->front_edge = or1->back_edge;
    or1->back_edge = e;
    return;
  }
  if( or1 )
  {
    if( e1 == or1->front_edge )
      or1->front_edge = e2;
    else
      or1->back_edge = e2;
  }
  if( or2 )
  {
    if( e2 == or2->front_edge )
      or2->front_edge = e1;
    else
      or2->back_edge = e1;
  }
  e1->outrec = or2;
  e2->outrec = or1;
  return;
}


/* FIXME TODO !!!PROBLEM!!! ~ We can NOT compute clockwiseness from area, fix me! */
static double clpGetOutPtArea( clpOutPt *op, int orientation_is_reversed )
{
  /* https://en.wikipedia.org/wiki/Shoelace_formula */
#if 0
  /* TODO FIXME ~ can't just call my polyarea stuff due to a linked list, need to expand here */
#else
  double area;
  clpOutPt *op2;
  area = 0.0;
  op2 = op;
  do
  {
    area += (double)( (op2->prev)->pt.y + op2->pt.y ) * (double)( (op2->prev)->pt.x - op2->pt.x );
    op2 = op2->next;
  } while( op2 != op );
  if( orientation_is_reversed )
    area = -area;
  return area * -0.5;
#endif
}

static inline double clpAreaTriangle( const clpPoint64 *pt1, const clpPoint64 *pt2, const clpPoint64 *pt3, int orientation_is_reversed )
{
  double area;
#if CLP_DEBUG_VERBOSE
  printf( "AreaTriangle %d %d ~ %d %d ~ %d %d\n", (int)pt1->x, (int)pt1->y, (int)pt2->x, (int)pt2->y, (int)pt3->x, (int)pt3->y );
#endif
  area = (double)( pt3->y + pt1->y ) * (double)( pt3->x - pt1->x ) +
    (double)( pt1->y + pt2->y ) * (double)( pt1->x - pt2->x ) +
    (double)( pt2->y + pt3->y ) * (double)( pt2->x - pt3->x );

  if( !( orientation_is_reversed ) )
    area = -area;
  return area;
}

static inline void clpReverseOutPts( clpOutPt *op )
{
  clpOutPt *op1, *op2;
  if( !op )
    return;
  op1 = op;
  do
  {
    op2 = op1->next;
    op1->next = op1->prev;
    op1->prev = op2;
    op1 = op2;
  } while( op1 != op );
  return;
}


static inline void clpSwapSides( clpOutRec *outrec )
{
  clpActiveEdge *e2;
  e2 = outrec->front_edge;
  outrec->front_edge = outrec->back_edge;
  outrec->back_edge = e2;
  outrec->pts = outrec->pts->next;
  return;
}


static inline clpOutRec *clpGetRealOutRec( clpOutRec *outrec )
{
  while( outrec && !outrec->pts )
    outrec = outrec->owner;
  return outrec;
}


static inline void clpUncoupleOutRec( clpActiveEdge *ae )
{
  clpOutRec *outrec;
  outrec = ae->outrec;
  if( !outrec )
    return;
  outrec->front_edge->outrec = 0;
  outrec->back_edge->outrec = 0;
  outrec->front_edge = 0;
  outrec->back_edge = 0;
  return;
}


static inline int clpAreReallyClose( const clpPoint64 *pt1, const clpPoint64 *pt2 )
{
#if 0
  return ( ( llabs( pt1->x - pt2->x ) < 2 ) && ( llabs( pt1->y - pt2->y ) < 2 ) );
#else
  return ( ( llabs( pt1->x - pt2->x ) | llabs( pt1->y - pt2->y ) ) < 2 );
#endif
}


static inline int clpIsValidClosedPath( const clpOutPt *op )
{
  /* also treat inconsequential polygons as invalid */
  return ( op && ( op->next != op ) && ( op->next != op->prev ) && ( ( op->next->next != op->prev ) || !( clpAreReallyClose( &op->pt, &(op->next)->pt ) && clpAreReallyClose( &op->pt, &(op->prev)->pt ) ) ) );
}


static inline int clpEdgesAdjacentInAEL( const clpIntersectNode *inode )
{
  return ( ( (inode->edge1)->aelnext == inode->edge2 ) || ( (inode->edge1)->aelprev == inode->edge2 ) );
}


////


static inline void clpDisposeAllOutRecs( clpClipper *clp )
{
  size_t index, count;
  clpOutRec *outrec;
  count = clpOutRecVector_GetCount( &clp->outreclist );
  for( index = 0 ; index < count ; index++ )
  {
    outrec = clpOutRecVector_Read( &clp->outreclist, index );
    if( outrec->pts )
      clpDisposeOutPts( clp, outrec );
    clpFreeOutRec( clp, outrec );
  }
  clpOutRecVector_Clear( &clp->outreclist );
  return;
}

static inline void clpDisposeVerticesAndLocalMinima( clpClipper *clp )
{
  size_t i, count;

  count = clpLocalMinimaVector_GetCount( &clp->minimalist );
  for( i = 0 ; i < count ; i++ )
    clpFreeLocalMinima( clp, clpLocalMinimaVector_Read( &clp->minimalist, i ) );
  clpLocalMinimaVector_Clear( &clp->minimalist );

  count = clpVertexList_GetCount( &clp->vertexlist );
  for( i = 0 ; i < count ; i++ )
    free( (void *)clpVertexList_Read( &clp->vertexlist, i ) );
  clpVertexList_Clear( &clp->vertexlist );

  return;
}


static void clpDeleteFromAEL( clpClipper *clp, clpActiveEdge *e )
{
  clpActiveEdge *prev, *next;
  prev = e->aelprev;
  next = e->aelnext;
#if 0
  if( !prev && !next && ( e != clp->activeedgelist ) )
    return;  /* Already deleted */
#endif
  if( prev )
    prev->aelnext = next;
  else
    clp->activeedgelist = next;
  if( next )
    next->aelprev = prev;
  clpFreeActiveEdge( clp, e );
  return;
}


static inline void clpCleanUp( clpClipper *clp )
{
  while( clp->activeedgelist )
    clpDeleteFromAEL( clp, clp->activeedgelist );
  mmHeap64sbFree( &clp->scanlinelist );
  clpIntersectNodeVector_Clear( &clp->intersectnodelist );
  clpDisposeAllOutRecs( clp );
  return;
}

static inline void clpClear( clpClipper *clp )
{
  clpCleanUp( clp );
  clpDisposeVerticesAndLocalMinima( clp );
  clp->minimalistwalkindex = 0;
  clp->minimalistsorted = 0;
  clp->openpathsflag = 0;
  return;
}


static inline void clpPushSortedScanline( clpClipper *clp, int64_t y )
{
#if CLP_DEBUG_VERBOSE
  printf( "&& InsertScanline %lld\n", (long long)y );
#endif
  mmHeap64sbSortedPush( &clp->scanlinelist, -y );
  return;
}

static inline void clpInsertScanline( clpClipper *clp, int64_t y )
{
#if CLP_DEBUG_VERBOSE
  printf( "&& InsertScanline %lld\n", (long long)y );
#endif
  mmHeap64sbInsert( &clp->scanlinelist, -y );
  return;
}


static inline int clpPopScanline( clpClipper *clp, int64_t *y )
{
  int64_t heapvalue, heaprep;
  if( !mmHeap64sbGetFirst( &clp->scanlinelist, &heapvalue ) )
    return 0;
  mmHeap64sbDeleteFirst( &clp->scanlinelist );
  for( ; ; )
  {
    if( !mmHeap64sbGetFirst( &clp->scanlinelist, &heaprep ) )
      break;
    if( heaprep != heapvalue )
      break;
    mmHeap64sbDeleteFirst( &clp->scanlinelist );
  }
  *y = -heapvalue;
#if CLP_DEBUG_VERBOSE
  printf( "&& Scanline GetFirst %llu\n", (long long)*y );
#endif
  return 1;
}


static inline void clpReset( clpClipper *clp )
{
  size_t mlindex, mlcount;
  int64_t pty, prevy;
  clpLocalMinima *lm;
  clpLocalMinima **lmlist;

  mlcount = clpLocalMinimaVector_GetCount( &clp->minimalist );
  lmlist = clpLocalMinimaVector_Access( &clp->minimalist, 0 );
  if( !clp->minimalistsorted )
  {
    /* Beyond a point, switch to radix sort */
#if CLP_CONFIG_ENABLE_RADIX_SORT
    if( mlcount < CLP_CONFIG_RADIX_SORT_COUNT_THRESHOLD && 0 )
      clpLocalMinimaHybridSort( lmlist, 0, mlcount, 128392130 );
    else
    {
      clpLocalMinimaRadixSortX( lmlist, 0, mlcount, 128 );
      clpLocalMinimaRadixSortY( lmlist, 0, mlcount, 128 );
    }
#else
    clpLocalMinimaHybridSort( lmlist, 0, mlcount, 128392130 );
#endif
    clp->minimalistsorted = 1;
  }

#if 0
  /* Inefficient insertion */
  for( mlindex = 0 ; mlindex < mlcount ; mlindex++ )
  {
    lm = clpLocalMinimaVector_Read( &clp->minimalist, mlindex );
    clpInsertScanline( clp, (lm->vertex)->pt.y );
  }
#else
  /* Our LocalMinima list is already sorted, push straight on the heap */
  prevy = INT64_MIN;
  for( mlindex = 0 ; mlindex < mlcount ; mlindex++ )
  {
    lm = clpLocalMinimaVector_Read( &clp->minimalist, mlindex );
    pty = (lm->vertex)->pt.y;
    if( pty != prevy )
    {
      clpPushSortedScanline( clp, pty );
      prevy = pty;
    }
  }
#endif

  clp->minimalistwalkindex = 0;
  clp->activeedgelist = 0;
  clp->sortededgelist = 0;
  clp->succeessflag = 1;
  return;
}


static void clpAddLocMin( clpClipper *clp, clpVertex *vertex, uint8_t polytype, uint8_t openflag )
{
  /* Make sure the vertex is added only once */
  if( vertex->flags & CLP_VERTEX_FLAG_LOCALMIN )
    return;
  vertex->flags |= CLP_VERTEX_FLAG_LOCALMIN;
  clpLocalMinimaVector_Push( &clp->minimalist, clpNewLocalMinima( clp, vertex, polytype, openflag ) );
  return;
}


static inline clpLocalMinima *clpPopLocalMinima( clpClipper *clp, int64_t y )
{
  clpLocalMinima *lm;
  if( clp->minimalistwalkindex == clpLocalMinimaVector_GetCount( &clp->minimalist ) )
    return 0;
  lm = clpLocalMinimaVector_Read( &clp->minimalist, clp->minimalistwalkindex );
  if( lm->vertex->pt.y != y )
    return 0;
  clp->minimalistwalkindex++;
  return lm;
}


static inline int clpIsContributingClosed( clpClipper *clp, const clpActiveEdge *e )
{
  switch( clp->fillrule )
  {
    case CLP_FILLRULE_EVENODD:
      break;
    case CLP_FILLRULE_NONZERO:
      if( abs( e->windingcount ) != 1 )
        return 0;
      break;
    default:
    {
      if( clp->fillrule == clp->fillpos )
      {
        if( e->windingcount != 1 )
          return 0;
      }
      else /* ( clp->fillrule == fillneg )*/
      {
        if( e->windingcount != -1 )
          return 0;
      }
    }
  }

  switch( clp->cliptype )
  {
    case CLP_CLIPTYPE_INTERSECTION:
      switch( clp->fillrule )
      {
        case CLP_FILLRULE_EVENODD:
        case CLP_FILLRULE_NONZERO:
          return ( e->windingcount2 != 0 );
      default:
        if( clp->fillrule == clp->fillpos )
          return ( e->windingcount2 > 0 );
        else
          return ( e->windingcount2 < 0 );
      }
      break;
    case CLP_CLIPTYPE_UNION:
      switch( clp->fillrule )
      {
        case CLP_FILLRULE_EVENODD:
        case CLP_FILLRULE_NONZERO:
          return ( e->windingcount2 == 0 );
        default:
          if( clp->fillrule == clp->fillpos )
            return ( e->windingcount2 <= 0 );
          else
            return ( e->windingcount2 >= 0 );
      }
      break;
    case CLP_CLIPTYPE_DIFFERENCE:
      if( clpGetPolyType( e ) == CLP_PATHTYPE_SUBJECT )
      {
        switch( clp->fillrule )
        {
          case CLP_FILLRULE_EVENODD:
          case CLP_FILLRULE_NONZERO:
            return ( e->windingcount2 == 0 );
          default:
            if( clp->fillrule == clp->fillpos )
              return ( e->windingcount2 <= 0 );
            else
              return ( e->windingcount2 >= 0 );
        }
      }
      else
      {
        switch( clp->fillrule )
        {
          case CLP_FILLRULE_EVENODD:
          case CLP_FILLRULE_NONZERO:
            return ( e->windingcount2 != 0 );
          default:
            if( clp->fillrule == clp->fillpos )
              return ( e->windingcount2 > 0 );
            else
              return ( e->windingcount2 < 0 );
        }
      }
      break;
    case CLP_CLIPTYPE_XOR:
      return 1;  /* XOR is always contributing unless open */
    default:
      return 0;
  }
  return 0;
}


static inline int clpIsContributingOpen( clpClipper *clp, const clpActiveEdge *e )
{
  int retval;
  switch( clp->cliptype )
  {
    case CLP_CLIPTYPE_INTERSECTION:
      retval = ( e->windingcount2 != 0 );
      break;
    case CLP_CLIPTYPE_UNION:
      retval = ( ( e->windingcount == 0 ) && ( e->windingcount2 == 0 ) );
      break;
    case CLP_CLIPTYPE_DIFFERENCE:
      retval = ( e->windingcount2 == 0 );
      break;
    case CLP_CLIPTYPE_XOR:
      retval = ( ( e->windingcount != 0 ) != ( e->windingcount2 != 0 ) );
      break;
    default:
      retval = 0;
      break;
  }
  return retval;
}


static inline void clpSetWindCountForClosedPathEdge( clpClipper *clp, clpActiveEdge *e )
{
  /* Wind counts refer to polygon regions not edges, so here an edge's WindCnt */
  /* indicates the higher of the wind counts for the two regions touching the */
  /* edge. (NB Adjacent regions can only ever have their wind counts differ by */
  /* one. Also, open paths have no meaningful wind directions or counts.) */
  clpActiveEdge *e2;
  uint8_t polytype;

  e2 = e->aelprev;
  /* Find the nearest closed path edge of the same PolyType in AEL (heading left) */
  polytype = clpGetPolyType( e );
  while( e2 && ( ( clpGetPolyType( e2 ) != polytype ) || clpIsActiveEdgeOpen( e2 ) ) )
    e2 = e2->aelprev;

  if( !e2 )
  {
    e->windingcount = e->windingdelta;
    e2 = clp->activeedgelist;
  }
  else if( clp->fillrule == CLP_FILLRULE_EVENODD )
  {
    e->windingcount = e->windingdelta;
    e->windingcount2 = e2->windingcount2;
    e2 = e2->aelnext;
  }
  else
  {
    /* NonZero, positive, or negative filling here */
    /* if e's windcnt is in the SAME direction as its WindDx, then polygon */
    /* filling will be on the right of 'e'. */
    /* NB neither e2.WindCnt nor e2.WindDx should ever be 0. */
    if( ( e2->windingcount * e2->windingdelta ) < 0 )
    {
      /* Opposite directions so 'e' is outside 'e2' */
      if( abs( e2->windingcount ) > 1 )
      {
        /* Outside prev poly but still inside another. */
        if( ( e2->windingdelta * e->windingdelta ) < 0 )
        {
          /* Reversing direction so use the same WC */
          e->windingcount = e2->windingcount;
        }
        else
        {
          /* Otherwise keep 'reducing' the WC by 1 (ie towards 0) */
          e->windingcount = e2->windingcount + e->windingdelta;
        }
      }
      else
      {
        /* Now outside all polys of same polytype so set own WC */
        e->windingcount = ( clpIsActiveEdgeOpen( e ) ? 1 : e->windingdelta );
      }
    }
    else
    {
      /* 'e' must be inside 'e2' */
      if( ( e2->windingdelta * e->windingdelta ) < 0 )
      {
        /* Reversing direction so use the same WC */
        e->windingcount = e2->windingcount;
      }
      else
      {
        /* Otherwise keep 'increasing' the WC by 1 (ie away from 0) */
        e->windingcount = e2->windingcount + e->windingdelta;
      }
    }
    e->windingcount2 = e2->windingcount2;
    e2 = e2->aelnext;  /* Get ready to calc WindCnt2 */
  }

  /* Update windingcount2 */
  if( clp->fillrule == CLP_FILLRULE_EVENODD )
  {
    while( e2 != e )
    {
      if( ( clpGetPolyType( e2 ) != polytype ) && !clpIsActiveEdgeOpen( e2 ) )
        e->windingcount2 = ( e->windingcount2 == 0 ? 1 : 0 );
      e2 = e2->aelnext;
    }
  }
  else
  {
    while( e2 != e )
    {
      if( ( clpGetPolyType( e2 ) != polytype ) && !clpIsActiveEdgeOpen( e2 ) )
        e->windingcount2 += e2->windingdelta;
      e2 = e2->aelnext;
    }
  }
  return;
}


static inline void clpSetWindCountForOpenPathEdge( clpClipper *clp, clpActiveEdge *e )
{
  int cnt1, cnt2;
  clpActiveEdge *e2;

  e2 = clp->activeedgelist;
  if( clp->fillrule == CLP_FILLRULE_EVENODD )
  {
    cnt1 = 0;
    cnt2 = 0;
    while( e2 != e )
    {
      if( clpGetPolyType( e2 ) == CLP_PATHTYPE_CLIP )
        cnt2++;
      else if( !clpIsActiveEdgeOpen( e2 ) )
        cnt1++;
      e2 = e2->aelnext;
    }
    e->windingcount = cnt1 & 1;
    e->windingcount2 = cnt2 & 1;
  }
  else
  {
    while( e2 != e )
    {
      if( clpGetPolyType( e2 ) == CLP_PATHTYPE_CLIP )
        e->windingcount2 += e2->windingdelta;
      else if( !clpIsActiveEdgeOpen( e2 ))
        e->windingcount += e2->windingdelta;
      e2 = e2->aelnext;
    }
  }
  return;
}


static int clpIsValidAelOrder( const clpActiveEdge *resident, const clpActiveEdge *newcomer )
{
  double d;
  int64_t y;
  int newcomerIsLeft;

  if( newcomer->scanbeamx != resident->scanbeamx )
    return ( newcomer->scanbeamx > resident->scanbeamx );

  /* Get the turning direction */
  d = clpCrossProduct( &resident->top, &newcomer->bot, &newcomer->top );
  if( d != 0.0 )
    return d < 0.0;

  /* Edges must be collinear to get here for starting open paths, place them according to the direction they're about to turn */
  if( !clpIsActiveEdgeMaxima( resident ) && ( resident->top.y > newcomer->top.y ) )
    return ( clpCrossProduct( &newcomer->bot, &resident->top, &clpGetNextVertex(resident)->pt ) <= 0.0 );
  else if( !clpIsActiveEdgeMaxima( newcomer ) && ( newcomer->top.y > resident->top.y ) )
    return ( clpCrossProduct( &newcomer->bot, &newcomer->top, &clpGetNextVertex(newcomer)->pt ) >= 0.0 );

  y = newcomer->bot.y;
  newcomerIsLeft = newcomer->leftboundflag;

  if( ( resident->bot.y != y ) || ( ((resident->local_min)->vertex)->pt.y != y ) )
    return newcomer->leftboundflag;
  else if( resident->leftboundflag != newcomerIsLeft ) /* Resident must also have just been inserted */
    return newcomerIsLeft;
  else if( clpCrossProduct( &clpGetPrevPrevVertex( resident )->pt, &resident->bot, &resident->top ) == 0.0 )
    return 1;
  else /* Compare turning direction of the alternate bound */
    return ( ( clpCrossProduct( &clpGetPrevPrevVertex( resident )->pt, &newcomer->bot, &clpGetPrevPrevVertex( newcomer )->pt ) > 0.0 ) == newcomerIsLeft );
}


static inline void clpInsertLeftEdge( clpClipper *clp, clpActiveEdge *e )
{
  clpActiveEdge *e2;
  if( !clp->activeedgelist )
  {
    e->aelprev = 0;
    e->aelnext = 0;
    clp->activeedgelist = e;
  }
  else if( !clpIsValidAelOrder( clp->activeedgelist, e ) )
  {
    e->aelprev = 0;
    e->aelnext = clp->activeedgelist;
    clp->activeedgelist->aelprev = e;
    clp->activeedgelist = e;
  }
  else
  {
    e2 = clp->activeedgelist;
    while( e2->aelnext && clpIsValidAelOrder( e2->aelnext, e ) )
      e2 = e2->aelnext;
    e->aelnext = e2->aelnext;
    if( e2->aelnext )
      e2->aelnext->aelprev = e;
    e->aelprev = e2;
    e2->aelnext = e;
  }
  return;
}


static inline void clpInsertRightEdge( clpActiveEdge *e, clpActiveEdge *e2 )
{
  e2->aelnext = e->aelnext;
  if( e->aelnext )
    e->aelnext->aelprev = e2;
  e2->aelprev = e;
  e->aelnext = e2;
  return;
}


static inline void clpPushHorz( clpClipper *clp, clpActiveEdge *e )
{
  e->selnext = clp->sortededgelist;
  clp->sortededgelist = e;
  return;
}


static int clpTestJoinWithPrev1( const clpActiveEdge *e, int64_t curr_y )
{
  /* This is marginally quicker than clpTestJoinWithPrev2 */
  /* but can only be used when e.PrevInAEL.currX is accurate */
  return clpIsHotEdge( e ) && !clpIsActiveEdgeOpen( e ) &&
    e->aelprev && ( (e->aelprev)->scanbeamx == e->scanbeamx ) &&
    clpIsHotEdge( e->aelprev ) && !clpIsActiveEdgeOpen( e->aelprev ) &&
    ( ( curr_y - e->top.y ) > 1 ) && ( ( curr_y - (e->aelprev)->top.y ) > 1 ) &&
    ( clpCrossProduct( &(e->aelprev)->top, &e->bot, &e->top ) == 0.0 );
}

static int clpTestJoinWithPrev2( const clpActiveEdge *e, const clpPoint64 *curr_pt )
{
  return clpIsHotEdge( e ) && !clpIsActiveEdgeOpen( e ) &&
    e->aelprev && !clpIsActiveEdgeOpen( e->aelprev ) &&
    clpIsHotEdge( e->aelprev ) && ( (e->aelprev)->top.y < e->bot.y ) &&
    ( llabs( clpGetActiveTopX( e->aelprev, curr_pt->y ) - curr_pt->x ) < 2 ) &&
    ( clpCrossProduct( &(e->aelprev)->top, curr_pt, &e->top ) == 0.0 );
}


static int clpTestJoinWithNext1( const clpActiveEdge *e, int64_t curr_y )
{
  /* This is marginally quicker than clpTestJoinWithNext2 */
  /* but can only be used when e.NextInAEL.currX is accurate */
  return clpIsHotEdge( e ) && !clpIsActiveEdgeOpen( e ) &&
    e->aelnext && ( (e->aelnext)->scanbeamx == e->scanbeamx ) &&
    clpIsHotEdge( e->aelnext ) && !clpIsActiveEdgeOpen( e->aelnext ) &&
    ( ( curr_y - e->top.y ) > 1 ) && ( ( curr_y - (e->aelnext)->top.y ) > 1 ) &&
    ( clpCrossProduct( &e->aelnext->top, &e->bot, &e->top ) == 0.0 );
}


static int clpTestJoinWithNext2( const clpActiveEdge *e, const clpPoint64 *curr_pt )
{
  return clpIsHotEdge( e ) && !clpIsActiveEdgeOpen( e ) &&
    e->aelnext && !clpIsActiveEdgeOpen( e->aelnext ) &&
    clpIsHotEdge( e->aelnext ) && ( (e->aelnext)->top.y < e->bot.y ) &&
    ( llabs( clpGetActiveTopX( e->aelnext, curr_pt->y ) - curr_pt->x ) < 2 ) &&
    ( clpCrossProduct( &(e->aelnext)->top, curr_pt, &e->top ) == 0.0 );
}


static inline int clpOutrecIsAscending( const clpActiveEdge *hotedge )
{
  return ( hotedge == (hotedge->outrec)->front_edge );
}


static clpOutPt *clpAddLocalMinPoly( clpClipper *clp, clpActiveEdge *e1, clpActiveEdge *e2, const clpPoint64 *pt, int is_new )
{
  clpOutRec *outrec;
  clpOutPt *op;
  clpActiveEdge *prevhotedge;

  outrec = clpNewOutRec( clp );
  outrec->idx = clpOutRecVector_GetCount( &clp->outreclist );
  clpOutRecVector_Push( &clp->outreclist, outrec );
  outrec->pts = 0;
  outrec->shapetree = 0;

  e1->outrec = outrec;
  e2->outrec = outrec;

  /* Setting the owner and inner/outer states (above) is an essential */
  /* precursor to setting edge 'sides' (ie left and right sides of output */
  /* polygons) and hence the orientation of output paths */

  if( clpIsActiveEdgeOpen( e1 ) )
  {
    /* WWW */
    outrec->owner = 0;
    outrec->openflag = 1;
    if( e1->windingdelta > 0 )
      clpSetSides( outrec, e1, e2 );
    else
      clpSetSides( outrec, e2, e1 );
  }
  else
  {
    prevhotedge = clpGetPrevHotEdge( e1 );
    /* e.windDx is the winding direction of the **input** paths and unrelated to the winding direction of output polygons */
    /* Output orientation is determined by e.outrec.frontE which is the ascending edge (see AddLocalMinPoly) */
    if( prevhotedge )
    {
      outrec->owner = prevhotedge->outrec;
      if( clpOutrecIsAscending( prevhotedge ) == is_new )
        clpSetSides( outrec, e2, e1 );
      else
        clpSetSides( outrec, e1, e2 );
    }
    else 
    {
      outrec->owner = 0;
      if( is_new )
        clpSetSides( outrec, e1, e2 );
      else
        clpSetSides( outrec, e2, e1 );
    }
  }
  op = clpNewOutPt( clp, pt, outrec );
  outrec->pts = op;
  return op;
}


static clpOutPt *clpAddOutPt( clpClipper *clp, const clpActiveEdge *e, const clpPoint64 pt )
{
  clpOutPt *new_op;
  clpOutRec *outrec;
  int to_front;
  clpOutPt *op_front;
  clpOutPt *op_back;
  /* outrec.outpts: a circular doubly-linked-list of outpt */
  outrec = e->outrec;
  to_front = clpIsFront( e );
  op_front = outrec->pts;
  op_back = op_front->next;
  if( to_front && clpPoint64Equal( &pt, &op_front->pt ) )
    new_op = op_front;
  else if( !to_front && clpPoint64Equal( &pt, &op_back->pt ) )
    new_op = op_back;
  else
  {
    new_op = clpNewOutPt( clp, &pt, outrec );
    op_back->prev = new_op;
    new_op->prev = op_front;
    new_op->next = op_back;
    op_front->next = new_op;
    if( to_front )
      outrec->pts = new_op;
  }
  return new_op;
}


static clpOutPt *clpStartOpenPath( clpClipper *clp, clpActiveEdge *e, const clpPoint64 *pt )
{
  clpOutPt *op;
  clpOutRec *outrec;
  outrec = clpNewOutRec( clp );
  outrec->idx = clpOutRecVector_GetCount( &clp->outreclist );
  clpOutRecVector_Push( &clp->outreclist, outrec );
  outrec->owner = 0;
  outrec->openflag = 1;
  outrec->pts = 0;
  outrec->shapetree = 0;

  if( e->windingdelta > 0 )
  {
    outrec->front_edge = e;
    outrec->back_edge = 0;
  }
  else
  {
    outrec->front_edge = 0;
    outrec->back_edge = e;
  }
  e->outrec = outrec;
  op = clpNewOutPt( clp, pt, outrec );
  outrec->pts = op;
  return op;
}


static inline clpActiveEdge *clpFindEdgeWithMatchingLocMin( clpActiveEdge *e )
{
  clpActiveEdge *result;
  result = e->aelnext;
  while( result )
  {
    if( result->local_min == e->local_min )
      return result;
    else if( !clpIsHorizontal( result ) && !clpPoint64Equal( &e->bot, &result->bot ) )
      result = 0;
    else
      result = result->aelnext;
  }
  result = e->aelprev;
  while( result )
  {
    if( result->local_min == e->local_min )
      return result;
    else if( !clpIsHorizontal( result ) && !clpPoint64Equal( &e->bot, &result->bot ) )
      return 0;
    else
      result = result->aelprev;
  }
  return result;
}


static inline clpJoiner *clpGetHorzTrialParent( clpClipper *clp, const clpOutPt *op )
{
  clpJoiner *joiner;
  joiner = op->joiner;
  while( joiner )
  {
    if( joiner->op1 == op )
    {
      if( joiner->next1 && ( (joiner->next1)->idx < 0 ) )
        return joiner;
      else
        joiner = joiner->next1;
    }
    else
    {
      if( joiner->next2 && ( (joiner->next2)->idx < 0 ) )
        return joiner;
      else
        joiner = joiner->next1;
    }
  }
  return joiner;
}


static inline int clpOutPtInTrialHorzList( clpClipper *clp, clpOutPt *op )
{
  return op->joiner && ( ( (op->joiner)->idx < 0 ) || clpGetHorzTrialParent( clp, op ) );
}


static inline clpJoiner *clpFindTrialJoinParent( clpJoiner **inoutjoiner, const clpOutPt *op )
{
  clpJoiner *joiner, *parent;
  joiner = *inoutjoiner;
  parent = joiner;
  while( parent )
  {
    if( op == parent->op1 )
    {
      if( parent->next1 && ( (parent->next1)->idx < 0 ) )
      {
        joiner = parent->next1;
        break;
      }
      parent = parent->next1;
    }
    else
    {
      if( parent->next2 && ( (parent->next2)->idx < 0 ) )
      {
        joiner = parent->next2;
        break;
      }
      parent = parent->next2;
    }
  }
  *inoutjoiner = joiner;
  return parent;
}


static void clpDeleteTrialHorzJoin( clpClipper *clp, clpOutPt *op )
{
  clpJoiner *joiner, *parenthorz, *parentop;

  if( !clp->horzjoinerlist )
    return;

  joiner = op->joiner;
  parentop = 0;
  while( joiner )
  {
    if( joiner->idx < 0 )
    {
      /* First remove joiner from FHorzTrials */
      if( joiner == clp->horzjoinerlist )
        clp->horzjoinerlist = joiner->nextH;
      else
      {
        parenthorz = clp->horzjoinerlist;
        while( parenthorz->nextH != joiner )
          parenthorz = parenthorz->nextH;
        parenthorz->nextH = joiner->nextH;
      }

      /* Now remove joiner from op's joiner list */
      if( !parentop )
      {
        /* Joiner must be first one in list */
        op->joiner = joiner->next1;
        clpFreeJoiner( clp, joiner );
        joiner = op->joiner;
      }
      else
      {
        /* The trial joiner isn't first */
        if( op == parentop->op1 )
          parentop->next1 = joiner->next1;
        else
          parentop->next2 = joiner->next1;
        clpFreeJoiner( clp, joiner );
        joiner = parentop;
      }
    }
    else
    {
      /* Not a trial join so look further along the linked list */
      parentop = clpFindTrialJoinParent( &joiner, op );
      if( !parentop )
        break;
    }
  }
  return;
}


static clpJoiner *clpFindJoinParent( const clpJoiner *joiner, clpOutPt *op )
{
  clpJoiner *result;
  result = op->joiner;
  for( ; ; )
  {
    if( op == result->op1 )
    {
      if( result->next1 == joiner )
        break;
      else
        result = result->next1;
    }
    else
    {
      if( result->next2 == joiner )
        break;
      else
        result = result->next2;
    }
  }
  return result;
}


/* This method deletes a single join, and it doesn't check for or delete trial horz joins */
static void clpDeleteJoin( clpClipper *clp, clpJoiner *joiner )
{
  clpOutPt *op1, *op2;
  clpJoiner *parentjoiner;

  op1 = joiner->op1;
  op2 = joiner->op2;

  if( op1->joiner != joiner )
  {
    parentjoiner = clpFindJoinParent( joiner, op1 );
    if( parentjoiner->op1 == op1 )
      parentjoiner->next1 = joiner->next1;
    else
      parentjoiner->next2 = joiner->next1;
  }
  else
    op1->joiner = joiner->next1;

  if( op2->joiner != joiner )
  {
    parentjoiner = clpFindJoinParent( joiner, op2 );
    if( parentjoiner->op1 == op2 )
      parentjoiner->next1 = joiner->next2;
    else
      parentjoiner->next2 = joiner->next2;
  }
  else
    op2->joiner = joiner->next2;

  clpJoinerVector_Write( &clp->joinerlist, joiner->idx, 0 );
  clpFreeJoiner( clp, joiner );
  return;
}

static void clpSafeDeleteOutPtJoiners( clpClipper *clp, clpOutPt* op )
{
  clpJoiner *joiner;

  for( joiner = op->joiner ; joiner ; joiner = op->joiner )
  {
    if( joiner->idx < 0 )
    {
      clpDeleteTrialHorzJoin( clp, op );
      continue;
    }
    else if( clp->horzjoinerlist )
    {
      if( clpOutPtInTrialHorzList( clp, joiner->op1 ) )
        clpDeleteTrialHorzJoin( clp, joiner->op1 );
      if( clpOutPtInTrialHorzList( clp, joiner->op2 ) )
        clpDeleteTrialHorzJoin( clp, joiner->op2 );
    }
    clpDeleteJoin( clp, joiner );
  }
  return;
}


static void clpSafeDisposeOutPts( clpClipper *clp, clpOutPt *op )
{
  clpOutRec *outrec;
  clpOutPt *tmp;
  outrec = clpGetRealOutRec( op->outrec );
  if( outrec->front_edge )
    (outrec->front_edge)->outrec = 0;
  if( outrec->back_edge )
    (outrec->back_edge)->outrec = 0;
  outrec->pts = 0;
  (op->prev)->next = 0;
  while( op )
  {
    clpSafeDeleteOutPtJoiners( clp, op );
    tmp = op->next;
    clpFreeOutPt( clp, op );
    op = tmp;
  }
  return;
}


static inline int clpValidateClosedPathEx( clpClipper *clp, clpOutRec *outrec )
{
  if( clpIsValidClosedPath( outrec->pts ) )
    return 1;
  clpSafeDisposeOutPts( clp, outrec->pts );
  return 0;
}


/* Check if two segments intersect, doesn't compute the actual intersection points */
#if 0
static inline int clpSegmentsIntersect( const clpPoint64 *seg1a, const clpPoint64 *seg1b, const clpPoint64 *seg2a, const clpPoint64 *seg2b )
{
  double dx1, dy1, dx2, dy2;
  double d2a1ax, d2a1ay, d2b1ax, d2b1ay, d2a1bx, d2a1by;
  double sma, smb, smc, smd;
  dx1 = (double)( seg1a->x - seg1b->x );
  dy1 = (double)( seg1a->y - seg1b->y );
  dx2 = (double)( seg2a->x - seg2b->x );
  dy2 = (double)( seg2a->y - seg2b->y );
  d2a1ax = (double)( seg2a->x - seg1a->x );
  d2a1ay = (double)( seg2a->y - seg1a->y );
  d2b1ax = (double)( seg2b->x - seg1a->x );
  d2b1ay = (double)( seg2b->y - seg1a->y );
  sma = ( ( dy1 * d2a1ax ) - ( dx1 * d2a1ay ) );
  smb = ( ( dy1 * d2b1ax ) - ( dx1 * d2b1ay ) );
  if( ( sma * smb ) >= 0.0 )
    return 0;
  d2a1bx = (double)( seg2a->x - seg1b->x );
  d2a1by = (double)( seg2a->y - seg1b->y );
  smc = ( ( dx2 * d2a1ay ) - ( dy2 * d2a1ax ) );
  smd = ( ( dx2 * d2a1by ) - ( dy2 * d2a1bx ) );
  if( ( smc * smd ) >= 0.0 )
    return 0;
  return 1;
/*
  return ( ( ( dy1 * ( seg2a->x - seg1a->x ) - dx1 * ( seg2a->y - seg1a->y ) ) *
    ( dy1 * ( seg2b->x - seg1a->x ) - dx1 * ( seg2b->y - seg1a->y ) ) < 0.0 ) &&
    ( ( dy2 * ( seg1a->x - seg2a->x ) - dx2 * ( seg1a->y - seg2a->y ) ) *
      ( dy2 * ( seg1b->x - seg2a->x ) - dx2 * ( seg1b->y - seg2a->y ) ) < 0.0 ) );
*/
}
#else
/* WWW */
static inline int clpSegmentsIntersect( const clpPoint64 *seg1a, const clpPoint64 *seg1b, const clpPoint64 *seg2a, const clpPoint64 *seg2b )
{
  return mathLineLine_Check( seg1a->x, seg1a->y, seg1b->x, seg1b->y, seg2a->x, seg2a->y, seg2b->x, seg2b->y );
}
#endif


static inline clpOutPt *clpDoSplitOp( clpClipper *clp, clpOutPt *outrecop, clpOutPt *splitop )
{
  clpOutPt *prevop;
  clpOutPt *nextnextop;
  clpOutPt *result;
  clpOutPt *newop, *newop2;
  clpOutRec *newoutrec;
  clpPointD ipd;
  clpPoint64 ip;
  double area1, area2, absarea2;

  prevop = splitop->prev;
  nextnextop = splitop->next->next;
  result = prevop;
  clpGetIntersectPoint( &prevop->pt, &splitop->pt, &splitop->next->pt, &nextnextop->pt, &ipd );
  ip = clpBuildPoint64( round( ipd.x ), round( ipd.y ) );
#if CLP_DEBUG_VERBOSE
  printf( "DoSplitOp\n" );
#endif
  area1 = clpGetOutPtArea( outrecop, clp->orientationreversedflag );
  area2 = clpAreaTriangle( &ip, &splitop->pt, &splitop->next->pt, clp->orientationreversedflag );
#if CLP_DEBUG_VERBOSE
  printf( "  Area %f %f\n", area1, area2 );
#endif
  if( clpPoint64Equal( &ip, &prevop->pt ) || clpPoint64Equal( &ip, &nextnextop->pt ) )
  {
    nextnextop->prev = prevop;
    prevop->next = nextnextop;
  }
  else
  {
    newop2 = clpNewOutPt( clp, &ip, prevop->outrec );
    newop2->prev = prevop;
    newop2->next = nextnextop;
    nextnextop->prev = newop2;
    prevop->next = newop2;
  }

  clpSafeDeleteOutPtJoiners( clp, splitop->next );
  clpSafeDeleteOutPtJoiners( clp, splitop );

  /* !!!PROBLEM!!! ~ We can NOT compute clockwiseness from area, fix me! */
  absarea2 = fabs( area2 );
  if( ( absarea2 >= 1.0 ) && ( ( ( absarea2 > fabs( area1 ) ) || ( ( area2 > 0 ) == ( area1 > 0 ) ) ) ) )
  {
    newoutrec = clpNewOutRec( clp );
    newoutrec->idx = clpOutRecVector_GetCount( &clp->outreclist );
    clpOutRecVector_Push( &clp->outreclist, newoutrec );
    newoutrec->owner = (prevop->outrec)->owner;
    newoutrec->shapetree = 0;
    splitop->outrec = newoutrec;
    splitop->next->outrec = newoutrec;

    newop = clpNewOutPt( clp, &ip, newoutrec );
    newop->prev = splitop->next;
    newop->next = splitop;
    newoutrec->pts = newop;
#if CLP_DEBUG_VERBOSE
    printf( "  DoSplitOp, %d %d\n", (int)newop->pt.x, (int)newop->pt.y );
#endif
    splitop->prev = newop;
    splitop->next->next = newop;
  }
  else
  {
    clpFreeOutPt( clp, splitop->next );
    clpFreeOutPt( clp, splitop );
  }
  return result;
}


static inline void clpFixSelfIntersects( clpClipper *clp, clpOutRec *outrec )
{
  clpOutPt *op2;
  op2 = outrec->pts;
  for( ; ; )
  {
    /* Three edged polygons can't self-intersect */
    if( op2->prev == (op2->next)->next )
      break;
    if( clpSegmentsIntersect( &(op2->prev)->pt, &op2->pt, &(op2->next)->pt, &((op2->next)->next)->pt ) )
    {
#if CLP_DEBUG_VERBOSE
      printf( "  FixSelfIntersects, %d %d\n", (int)op2->pt.x, (int)op2->pt.y );
#endif
      if( ( op2 == outrec->pts ) || ( op2->next == outrec->pts ) )
        outrec->pts = (outrec->pts)->prev;
      op2 = clpDoSplitOp( clp, outrec->pts, op2 );
      outrec->pts = op2;
      continue;
    }
    else
      op2 = op2->next;
    if( op2 == outrec->pts )
      break;
  }
  return;
}


static void clpCleanCollinear( clpClipper *clp, clpOutRec *outrec )
{
  clpOutPt *startOp, *op2;

  outrec = clpGetRealOutRec( outrec );
  if( !outrec || outrec->openflag || outrec->front_edge || !clpValidateClosedPathEx( clp, outrec ) )
    return;

  startOp = outrec->pts;

#if CLP_DEBUG_VERBOSE
  DebugDebugPrintPtChain( startOp, "CleanCollinear" );
#endif

  op2 = startOp;
  for( ; ; )
  {
    if( op2->joiner )
      return;
#if CLP_DEBUG_VERBOSE
    printf( "  CleanCollinear, %d %d\n", (int)op2->pt.x, (int)op2->pt.y );
#endif
    /* NB if preservecollinear == true, then only remove 180 deg. spikes */
    if( ( clpCrossProduct( &op2->prev->pt, &op2->pt, &op2->next->pt ) == 0.0 ) && ( clpPoint64Equal( &op2->pt, &op2->prev->pt ) || clpPoint64Equal( &op2->pt, &op2->next->pt ) || !clp->preservecollinearflag || clpDotProduct( &(op2->prev)->pt, &op2->pt, &(op2->next)->pt ) < 0 ) )
    {
      if( op2 == outrec->pts )
        outrec->pts = op2->prev;
      op2 = clpDisposeOutPt( clp, op2 );
      if( !clpValidateClosedPathEx( clp, outrec ) )
      {
        outrec->pts = 0;
        return;
      }
      startOp = op2;
      continue;
    }
    op2 = op2->next;
    if( op2 == startOp )
      break;
  }
  clpFixSelfIntersects( clp, outrec );
  return;
}


static void clpJoinOutrecPaths( clpActiveEdge *e1, clpActiveEdge *e2 )
{
  /* Join e2 outrec path onto e1 outrec path and then delete e2 outrec path */
  /* pointers (only very rarely do the joining ends share the same coords) */
  clpOutPt *p1_st, *p2_st, *p1_end, *p2_end;

  p1_st = (e1->outrec)->pts;
  p2_st = (e2->outrec)->pts;
  p1_end = p1_st->next;
  p2_end = p2_st->next;
  if( clpIsFront( e1 ) )
  {
    p2_end->prev = p1_st;
    p1_st->next = p2_end;
    p2_st->next = p1_end;
    p1_end->prev = p2_st;
    e1->outrec->pts = p2_st;
    e1->outrec->front_edge = (e2->outrec)->front_edge;
    if( (e1->outrec)->front_edge )
      ((e1->outrec)->front_edge)->outrec = e1->outrec;
  }
  else
  {
    p1_end->prev = p2_st;
    p2_st->next = p1_end;
    p1_st->next = p2_end;
    p2_end->prev = p1_st;
    (e1->outrec)->back_edge = (e2->outrec)->back_edge;
    if( (e1->outrec)->back_edge )
      ((e1->outrec)->back_edge)->outrec = e1->outrec;
  }

  /* An owner must have a lower idx, otherwise it can't be a valid owner */
  if( (e2->outrec)->owner && ( ((e2->outrec)->owner)->idx < (e1->outrec)->idx ) )
  {
    if( !(e1->outrec)->owner || ( ((e2->outrec)->owner)->idx < ((e1->outrec)->owner)->idx) )
      (e1->outrec)->owner = (e2->outrec)->owner;
  }

  /* After joining, the e2->outrec must contains no vertices */
  (e2->outrec)->front_edge = 0;
  (e2->outrec)->back_edge = 0;
  (e2->outrec)->pts = 0;
  (e2->outrec)->owner = e1->outrec;

  if( clpIsActiveEdgeOpenEnd( e1 ) )
  {
    (e2->outrec)->pts = (e1->outrec)->pts;
    (e1->outrec)->pts = 0;
  }

  /* Both e1 and e2 are maxima and are about to be dropped from the Actives list */
  e1->outrec = 0;
  e2->outrec = 0;
  return;
}


static clpOutPt *clpAddLocalMaxPoly( clpClipper *clp, clpActiveEdge *e1, clpActiveEdge *e2, const clpPoint64 pt )
{
  clpOutPt *result;
  clpOutRec *outrec;
  if( clpIsFront( e1 ) == clpIsFront( e2 ) )
  {
    if( clpIsActiveEdgeOpen( e1 ) )
    {
      if( e1->windingdelta < 0 )
        clpSwapSides( e2->outrec );
      else
        clpSwapSides( e1->outrec );
    }
    else
    {
      clp->succeessflag = 0;
      return 0;
    }
  }

  result = clpAddOutPt( clp, e1, pt );
  if( e1->outrec == e2->outrec )
  {
    outrec = e1->outrec;
    outrec->pts = result;
    clpUncoupleOutRec( e1 );
    if( !clpIsActiveEdgeOpen( e1 ) )
      clpCleanCollinear( clp, outrec );
    result = outrec->pts;
    if( clp->polytreeflag && outrec->owner && !(outrec->owner)->front_edge)
      outrec->owner = clpGetRealOutRec( (outrec->owner)->owner );
  }
  /* and to preserve the winding orientation of outrec */
  else if( clpIsActiveEdgeOpen( e1 ) )
  {
    if( e1->windingdelta < 0 )
      clpJoinOutrecPaths( e1, e2 );
    else
      clpJoinOutrecPaths( e2, e1 );
  }
  else if( e1->outrec->idx < e2->outrec->idx )
    clpJoinOutrecPaths( e1, e2 );
  else
    clpJoinOutrecPaths( e2, e1 );
  return result;
}


static void clpAddJoin( clpClipper *clp, clpOutPt *op1, clpOutPt *op2 )
{
  clpJoiner *j;
  /* Unless op1->next or op1->prev crosses the start-end divide, don't waste time trying to join adjacent vertices */
  if( ( op1->outrec == op2->outrec ) && ( ( op1 == op2 ) || ( ( op1->next == op2 ) && ( op1 != op1->outrec->pts ) ) || ( ( op2->next == op1 ) && ( op2 != op1->outrec->pts ) ) ) )
    return;
  j = clpNewJoiner( clp, op1, op2, 0 );
  j->idx = clpJoinerVector_GetCount( &clp->joinerlist );
#if CLP_DEBUG_VERBOSE
  printf( "  AddJoin %d ; %d,%d %d,%d\n", j->idx, (int)op1->pt.x, (int)op1->pt.y, (int)op2->pt.x, (int)op2->pt.y );
#endif
  clpJoinerVector_Push( &clp->joinerlist, j );
  return;
}


static inline int clpIsZeroOrOne( int count )
{
#if 0
  return ( ( count == 0 ) || ( count == 1 ) );
#elif 0
  return ( ( count & ~1 ) == 0 );
#else
  return ( (unsigned int)count < 2 );
#endif
}

static clpOutPt *clpIntersectEdges( clpClipper *clp, clpActiveEdge *e1, clpActiveEdge *e2, const clpPoint64 pt )
{
  int swapwc, e1wc, e2wc;
  int e1_wc0or1, e2_wc0or1;
  clpActiveEdge *edge_o, *edge_c;
  clpOutPt *newoutpt;
  clpActiveEdge *e3;
  clpOutPt *outptlm;
  int64_t e1wc2, e2wc2;
  clpLocalMinima *e1lm, *e2lm;

  /* First handle open paths, if any ~ could be put in separate function? */
  if( clp->openpathsflag && ( clpIsActiveEdgeOpen( e1 ) || clpIsActiveEdgeOpen( e2 ) ) )
  {
    if( clpIsActiveEdgeOpen( e1 ) && clpIsActiveEdgeOpen( e2 ) )
      return 0;

    if( clpIsActiveEdgeOpen( e1 ) )
    {
      edge_o = e1;
      edge_c = e2;
    }
    else
    {
      edge_o = e2;
      edge_c = e1;
    }

    if( abs( edge_c->windingcount ) != 1 )
      return 0;

    if( clp->cliptype == CLP_CLIPTYPE_UNION )
    {
      if( ( abs( e2->windingcount ) != 1 ) || ( e2->windingcount2 != 0 ) )
        return 0;
    }
    else if( clpIsSamePolyType( e1, e2 ) || ( abs( e2->windingcount ) != 1 ) )
      return 0;

    /* Toggle contribution */
    if( clpIsHotEdge( edge_o ) )
    {
      newoutpt = clpAddOutPt( clp, edge_o, pt );
      edge_o->outrec = 0;
      return newoutpt;
    }
    else if( clpPoint64Equal( &pt, &(edge_o->local_min)->vertex->pt ) && !clpIsVertexOpenEnd( (edge_o->local_min)->vertex ) ) // horizontal edges can pass under open paths at a local minima
    {
      /* Find the other side of the localminima and if it's 'hot' join up with it */
      e3 = clpFindEdgeWithMatchingLocMin( edge_o );
      if( clpIsHotEdge( e3 ) )
      {
        edge_o->outrec = e3->outrec;
        if( edge_o->windingdelta > 0 )
          clpSetSides( e3->outrec, edge_o, e3 );
        else
          clpSetSides( e3->outrec, e3, edge_o );
        return e3->outrec->pts;
      }
      else
        return clpStartOpenPath( clp, edge_o, &pt );
    }
    else
      return clpStartOpenPath( clp, edge_o, &pt );
  }

  e1lm = e1->local_min;
  e2lm = e2->local_min;

  /* Closed path handling, update winding counts */
  if( ( e1lm->polytype == e2lm->polytype ) )
  {
    if( clp->fillrule == CLP_FILLRULE_EVENODD )
    {
      swapwc = e1->windingcount;
      e1->windingcount = e2->windingcount;
      e2->windingcount = swapwc;
    }
    else
    {
      if( ( e1->windingcount + e2->windingdelta ) == 0 )
        e1->windingcount = -e1->windingcount;
      else
        e1->windingcount += e2->windingdelta;
      if( ( e2->windingcount - e1->windingdelta ) == 0 )
        e2->windingcount = -e2->windingcount;
      else
        e2->windingcount -= e1->windingdelta;
    }
  }
  else
  {
    if( clp->fillrule != CLP_FILLRULE_EVENODD )
    {
      e1->windingcount2 += e2->windingdelta;
      e2->windingcount2 -= e1->windingdelta;
    }
    else
    {
      e1->windingcount2 = ( e1->windingcount2 == 0 ? 1 : 0 );
      e2->windingcount2 = ( e2->windingcount2 == 0 ? 1 : 0 );
    }
  }

  switch( clp->fillrule )
  {
    case CLP_FILLRULE_EVENODD:
    case CLP_FILLRULE_NONZERO:
      e1wc = abs( e1->windingcount );
      e2wc = abs( e2->windingcount );
      break;
    default:
      if( clp->fillrule == clp->fillpos )
      {
        e1wc = e1->windingcount;
        e2wc = e2->windingcount;
      }
      else
      {
        e1wc = -e1->windingcount;
        e2wc = -e2->windingcount;
      }
      break;
  }

  e1_wc0or1 = clpIsZeroOrOne( e1wc );
  e2_wc0or1 = clpIsZeroOrOne( e2wc );

  if( ( !clpIsHotEdge( e1 ) && !e1_wc0or1 ) || ( !clpIsHotEdge( e2 ) && !e2_wc0or1 ) )
    return 0;

  /* Process the intersection */
  newoutpt = 0;
  if( clpIsHotEdge( e1 ) && clpIsHotEdge( e2 ) )
  {
    if( !e1_wc0or1 || !e2_wc0or1 || ( ( e1lm->polytype != e2lm->polytype ) && ( clp->cliptype != CLP_CLIPTYPE_XOR ) ) )
      newoutpt = clpAddLocalMaxPoly( clp, e1, e2, pt );
    else if( clpIsFront( e1 ) || ( e1->outrec == e2->outrec ) )
    {
      newoutpt = clpAddLocalMaxPoly( clp, e1, e2, pt );
      outptlm = clpAddLocalMinPoly( clp, e1, e2, &pt, 0 );
      if( newoutpt && clpPoint64Equal( &newoutpt->pt, &outptlm->pt ) && !clpIsHorizontal( e1 ) && !clpIsHorizontal( e2 ) && ( clpCrossProduct( &e1->bot, &newoutpt->pt, &e2->bot ) == 0.0 ) )
        clpAddJoin( clp, newoutpt, outptlm );
    }
    else
    {
      newoutpt = clpAddOutPt( clp, e1, pt );
      clpAddOutPt( clp, e2, pt );
      clpSwapOutrecs( e1, e2 );
    }
  }
  else if( clpIsHotEdge( e1 ) )
  {
    newoutpt = clpAddOutPt( clp, e1, pt );
    clpSwapOutrecs( e1, e2 );
  }
  else if( clpIsHotEdge( e2 ) )
  {
    newoutpt = clpAddOutPt( clp, e2, pt );
    clpSwapOutrecs( e1, e2 );
  }
  else
  {
    switch( clp->fillrule )
    {
      case CLP_FILLRULE_EVENODD:
      case CLP_FILLRULE_NONZERO:
        e1wc2 = abs( e1->windingcount2 );
        e2wc2 = abs( e2->windingcount2 );
        break;
      default:
        if( clp->fillrule == clp->fillpos )
        {
          e1wc2 = e1->windingcount2;
          e2wc2 = e2->windingcount2;
        }
        else
        {
          e1wc2 = -e1->windingcount2;
          e2wc2 = -e2->windingcount2;
        }
        break;
    }
    if( e1lm->polytype != e2lm->polytype )
      newoutpt = clpAddLocalMinPoly( clp, e1, e2, &pt, 0 );
    else if( ( e1wc == 1 ) && ( e2wc == 1 ) )
    {
      newoutpt = 0;
      switch( clp->cliptype )
      {
        case CLP_CLIPTYPE_UNION:
          if( ( e1wc2 <= 0 ) && ( e2wc2 <= 0 ) )
            newoutpt = clpAddLocalMinPoly( clp, e1, e2, &pt, 0 );
          break;
        case CLP_CLIPTYPE_DIFFERENCE:
          if( ( ( e1lm->polytype == CLP_PATHTYPE_CLIP ) && ( e1wc2 > 0 ) && ( e2wc2 > 0 ) ) || ( ( e1lm->polytype == CLP_PATHTYPE_SUBJECT ) && ( e1wc2 <= 0 ) && ( e2wc2 <= 0 ) ) )
            newoutpt = clpAddLocalMinPoly( clp, e1, e2, &pt, 0 );
          break;
        case CLP_CLIPTYPE_XOR:
          newoutpt = clpAddLocalMinPoly( clp, e1, e2, &pt, 0 );
          break;
        default:
          if( ( e1wc2 > 0 ) && ( e2wc2 > 0 ) )
            newoutpt = clpAddLocalMinPoly( clp, e1, e2, &pt, 0 );
          break;
      }
    }
  }
  return newoutpt;
}


static void clpSwapPositionsInAEL( clpClipper *clp, clpActiveEdge *e1, clpActiveEdge *e2 )
{
  clpActiveEdge *next, *prev;
  /* preconditon: e1 must be immediately to the left of e2 */
  next = e2->aelnext;
  if( next )
    next->aelprev = e1;
  prev = e1->aelprev;
  if( prev )
    prev->aelnext = e2;
  e2->aelprev = prev;
  e2->aelnext = e1;
  e1->aelprev = e2;
  e1->aelnext = next;
  if( !e2->aelprev )
    clp->activeedgelist = e2;
  return;
}


static inline void clpInsertLocalMinimaIntoAEL( clpClipper *clp, int64_t scanbeamboty )
{
  clpLocalMinima *localminima;
  clpActiveEdge *leftbound, *rightbound;
  int contributing;
  clpOutPt *op;

  /* Add any local minima (if any) at scambeamboty */
  /* Horizontal local minima edges should contain (lm->vertex)->prev */
  while( ( localminima = clpPopLocalMinima( clp, scanbeamboty ) ) )
  {
    if( (localminima->vertex)->flags & CLP_VERTEX_FLAG_OPENSTART )
      leftbound = 0;
    else
    {
      leftbound = clpNewActiveEdge( clp );
      leftbound->bot = (localminima->vertex)->pt;
      leftbound->scanbeamx = leftbound->bot.x;
      leftbound->windingdelta = -1,
      leftbound->vertextop = (localminima->vertex)->prev;  //ie descending
      leftbound->top = (leftbound->vertextop)->pt;
      leftbound->outrec = 0;
      leftbound->local_min = localminima;
      clpActiveUpdateDx( leftbound );
    }

    if( (localminima->vertex)->flags & CLP_VERTEX_FLAG_OPENEND )
      rightbound = 0;
    else
    {
      rightbound = clpNewActiveEdge( clp );
      rightbound->bot = (localminima->vertex)->pt;
      rightbound->scanbeamx = rightbound->bot.x;
      rightbound->windingdelta = 1,
      rightbound->vertextop = (localminima->vertex)->next;  //ie ascending
      rightbound->top = (rightbound->vertextop)->pt;
      rightbound->outrec = 0;
      rightbound->local_min = localminima;
      clpActiveUpdateDx( rightbound );
    }

    /* Currently LeftB is just the descending bound and RightB is the ascending. */
    /* Now if the LeftB isn't on the left of RightB then we need swap them. */
    if( leftbound && rightbound )
    {
      if( clpIsHorizontal( leftbound ) )
      {
        if( clpIsHeadingRightHorz( leftbound ) )
          clpSwapActives( &leftbound, &rightbound );
      }
      else if( clpIsHorizontal( rightbound ) )
      {
        if( clpIsHeadingLeftHorz( rightbound ) )
          clpSwapActives( &leftbound, &rightbound );
      }
      else if( leftbound->slopedx < rightbound->slopedx )
        clpSwapActives( &leftbound, &rightbound );
    }
    else if( !leftbound )
    {
      leftbound = rightbound;
      rightbound = 0;
    }

    leftbound->leftboundflag = 1;
    clpInsertLeftEdge( clp, leftbound );

    if( clpIsActiveEdgeOpen( leftbound ) )
    {
      clpSetWindCountForOpenPathEdge( clp, leftbound );
      contributing = clpIsContributingOpen( clp, leftbound );
    }
    else
    {
      clpSetWindCountForClosedPathEdge( clp, leftbound );
      contributing = clpIsContributingClosed( clp, leftbound );
    }

    if( rightbound )
    {
      rightbound->leftboundflag = 0;
      rightbound->windingcount = leftbound->windingcount;
      rightbound->windingcount2 = leftbound->windingcount2;
      clpInsertRightEdge( leftbound, rightbound );
      if( contributing )
      {
        clpAddLocalMinPoly( clp, leftbound, rightbound, &leftbound->bot, 1 );
        if( !clpIsHorizontal( leftbound ) && clpTestJoinWithPrev1( leftbound, scanbeamboty ) )
        {
          op = clpAddOutPt( clp, leftbound->aelprev, leftbound->bot );
          clpAddJoin( clp, op, (leftbound->outrec)->pts );
        }
      }

      while( rightbound->aelnext && clpIsValidAelOrder( rightbound->aelnext, rightbound ) )
      {
        clpIntersectEdges( clp, rightbound, rightbound->aelnext, rightbound->bot );
        clpSwapPositionsInAEL( clp, rightbound, rightbound->aelnext );
      }

      if( !clpIsHorizontal( rightbound ) && clpTestJoinWithNext1( rightbound, scanbeamboty ) )
      {
        op = clpAddOutPt( clp, rightbound->aelnext, rightbound->bot );
        clpAddJoin( clp, rightbound->outrec->pts, op );
      }

      if( clpIsHorizontal( rightbound ) )
        clpPushHorz( clp, rightbound );
      else
        clpInsertScanline( clp, rightbound->top.y );
    }
    else if( contributing )
      clpStartOpenPath( clp, leftbound, &leftbound->bot );

    if( clpIsHorizontal( leftbound ) )
      clpPushHorz( clp, leftbound );
    else
      clpInsertScanline( clp, leftbound->top.y );
  }
  return;
}


static inline clpActiveEdge *clpPopHorz( clpClipper *clp )
{
  clpActiveEdge *ae;
  ae = clp->sortededgelist;
  if( !ae )
    return 0;
  clp->sortededgelist = (clp->sortededgelist)->selnext;
  return ae;
}


static inline void clpUpdateOutrecOwner( clpOutRec *outrec )
{
  clpOutPt *opCurr;
  opCurr = outrec->pts;
  for( ; ; )
  {
    opCurr->outrec = outrec;
    opCurr = opCurr->next;
    if( opCurr == outrec->pts )
      return;
  }
  return;
}


static void clpCompleteSplit( clpClipper *clp, clpOutPt *op1, clpOutPt *op2, clpOutRec *outrec )
{
  double area1, area2;
  int signs_change;
  area1 = clpGetOutPtArea( op1, clp->orientationreversedflag );
  area2 = clpGetOutPtArea( op2, clp->orientationreversedflag );
  signs_change = ( ( area1 > 0.0 ) == ( area2 < 0.0 ) );

  /* FIXME TODO !!!PROBLEM!!! ~ Area math is unreliable */

  if( ( area1 == 0.0 ) || ( signs_change && ( fabs( area1 ) < 2.0 ) ) )
  {
    clpSafeDisposeOutPts( clp, op1 );
    outrec->pts = op2;
  }
  else if( ( area2 == 0.0 ) || ( signs_change && ( fabs( area2 ) < 2.0 ) ) )
  {
    clpSafeDisposeOutPts( clp, op2 );
    outrec->pts = op1;
  }
  else
  {
    clpOutRec *newOr = clpNewOutRec( clp );
    newOr->idx = clpOutRecVector_GetCount( &clp->outreclist );
    clpOutRecVector_Push( &clp->outreclist, newOr );
    newOr->shapetree = 0;

    if( clp->polytreeflag )
    {
#if CLP_DEBUG_VERBOSE
      printf( "## PushSplit : Outrec.idx %d\n", (int)outrec->idx );
      if( newOr->pts )
        DebugDebugPrintPtChain( newOr->pts, "   OutRecSplit" );
#endif
      clpOutRecVector_Push( &outrec->splits, newOr );
    }

    if( fabs( area1 ) >= fabs( area2 ) )
    {
      outrec->pts = op1;
      newOr->pts = op2;
    }
    else
    {
      outrec->pts = op2;
      newOr->pts = op1;
    }

    /* FIXME TODO !!!PROBLEM!!! ~ We can NOT compute clockwiseness from area, fix me! */
    if( ( area1 > 0.0 ) == ( area2 > 0.0 ) )
      newOr->owner = outrec->owner;
    else
      newOr->owner = outrec;

    clpUpdateOutrecOwner( newOr );
    clpCleanCollinear( clp, newOr );
  }
  return;
}


static void clpUpdateEdgeIntoAEL( clpClipper *clp, clpActiveEdge *e )
{
  clpOutPt *op1, *op2;
  e->bot = e->top;
  e->vertextop = clpGetNextVertex( e );
  e->top = e->vertextop->pt;
  e->scanbeamx = e->bot.x;
  clpActiveUpdateDx( e );
  if( clpIsHorizontal( e ) )
    return;
  clpInsertScanline( clp, e->top.y );
  if( clpTestJoinWithPrev1( e, e->bot.y ) )
  {
    op1 = clpAddOutPt( clp, e->aelprev, e->bot );
    op2 = clpAddOutPt( clp, e, e->bot );
    clpAddJoin( clp, op1, op2 );
  }
  return;
}


static inline void clpAdjustCurrXAndCopyToSEL( clpClipper *clp, const int64_t scantopy )
{
  clpActiveEdge *e;
  e = clp->activeedgelist;
  clp->sortededgelist = e;
  while( e )
  {
    e->selprev = e->aelprev;
    e->selnext = e->aelnext;
    e->jump = e->selnext;
    e->scanbeamx = clpGetActiveTopX( e, scantopy );
    e = e->aelnext;
  }
  return;
}


#if 0
static inline int clpVertexLineLineIntersection( int64_t *hitpt, int64_t l0p0x, int64_t l0p0y, int64_t l0p1x, int64_t l0p1y, int64_t l1p0x, int64_t l1p0y, int64_t l1p1x, int64_t l1p1y )
{
  int64_t ln0a, ln0b, ln1a, ln1b;
  double ln0c, ln1c, det;
  ln0a = l0p1y - l0p0y;
  ln0b = l0p0x - l0p1x;
  ln1a = l1p1y - l1p0y;
  ln1b = l1p0x - l1p1x;
  ln0c = -( (double)ln0a * (double)l0p0x ) - ( (double)ln0b * (double)l0p0y );
  ln1c = -( (double)ln1a * (double)l1p0x ) - ( (double)ln1b * (double)l1p0y );
  det = ( (double)ln1a * (double)ln0b ) - ( (double)ln0a * (double)ln1b );
  if( det != 0.0 )
  {
    hitpt[0] = (int64_t)nearbyint( ( ( (double)ln1b * ln0c ) - ( (double)ln0b * ln1c ) ) / det );
    hitpt[1] = (int64_t)nearbyint( ( ( (double)ln0a * ln1c ) - ( (double)ln1a * ln0c ) ) / det );
    return 1;
  }
  return 0;
}
#elif 0
static inline int clpVertexLineLineIntersection( int64_t *hitpt, int64_t l0p0x, int64_t l0p0y, int64_t l0p1x, int64_t l0p1y, int64_t l1p0x, int64_t l1p0y, int64_t l1p1x, int64_t l1p1y )
{
  int64_t ln0a, ln0b, ln1a, ln1b;
  int64_t originx, originy;
  double ln0c, ln1c, det;
  ln0a = l0p1y - l0p0y;
  ln0b = l0p0x - l0p1x;
  ln1a = l1p1y - l1p0y;
  ln1b = l1p0x - l1p1x;
  originx = (int64_t)( 0.25 * ( (double)l0p0x + (double)l0p1x + (double)l1p0x + (double)l1p1x ) );
  originy = (int64_t)( 0.25 * ( (double)l0p0y + (double)l0p1y + (double)l1p0y + (double)l1p1y ) );
  ln0c = -( (double)ln0a * (double)( l0p0x - originx ) ) - ( (double)ln0b * (double)( l0p0y - originy ) );
  ln1c = -( (double)ln1a * (double)( l1p0x - originx ) ) - ( (double)ln1b * (double)( l1p0y - originy ) );
  det = ( (double)ln1a * (double)ln0b ) - ( (double)ln0a * (double)ln1b );
  if( det != 0.0 )
  {
    hitpt[0] = originx + (int64_t)nearbyint( ( ( (double)ln1b * ln0c ) - ( (double)ln0b * ln1c ) ) / det );
    hitpt[1] = originy + (int64_t)nearbyint( ( ( (double)ln0a * ln1c ) - ( (double)ln1a * ln0c ) ) / det );
    return 1;
  }
  return 0;
}
#elif 0
static inline int clpVertexLineLineIntersection_C( int64_t *hitpt, int64_t l0p0x, int64_t l0p0y, int64_t l0p1x, int64_t l0p1y, int64_t l1p0x, int64_t l1p0y, int64_t l1p1x, int64_t l1p1y )
{
  int64_t ln0a, ln0b, ln1a, ln1b;
  int64_t originx, originy;
  double ln0c, ln1c, det;
  ln0a = l0p1y - l0p0y;
  ln0b = l0p0x - l0p1x;
  ln1a = l1p1y - l1p0y;
  ln1b = l1p0x - l1p1x;
  det = ( (double)ln1a * (double)ln0b ) - ( (double)ln0a * (double)ln1b );
  if( det == 0.0 )
    return 0;
  ln0c = -( (double)ln0a * (double)l0p0x ) - ( (double)ln0b * (double)l0p0y );
  ln1c = -( (double)ln1a * (double)l1p0x ) - ( (double)ln1b * (double)l1p0y );
  originx = (int64_t)nearbyint( ( ( (double)ln1b * ln0c ) - ( (double)ln0b * ln1c ) ) / det );
  originy = (int64_t)nearbyint( ( ( (double)ln0a * ln1c ) - ( (double)ln1a * ln0c ) ) / det );
  det = ( (double)ln1a * (double)ln0b ) - ( (double)ln0a * (double)ln1b );
  if( det != 0.0 )
  {
    ln0c = -( (double)ln0a * (double)( l0p0x - originx ) ) - ( (double)ln0b * (double)( l0p0y - originy ) );
    ln1c = -( (double)ln1a * (double)( l1p0x - originx ) ) - ( (double)ln1b * (double)( l1p0y - originy ) );
    hitpt[0] = originx + (int64_t)nearbyint( ( ( (double)ln1b * ln0c ) - ( (double)ln0b * ln1c ) ) / det );
    hitpt[1] = originy + (int64_t)nearbyint( ( ( (double)ln0a * ln1c ) - ( (double)ln1a * ln0c ) ) / det );
    return 1;
  }
  return 0;
}
#else
/* WWW */
static inline int clpVertexLineLineIntersection_C( int64_t *hitpt, int64_t l0p0x, int64_t l0p0y, int64_t l0p1x, int64_t l0p1y, int64_t l1p0x, int64_t l1p0y, int64_t l1p1x, int64_t l1p1y )
{
  return mathLineLine_HitAnySafeRobust( hitpt, l0p0x, l0p0y, l0p1x, l0p1y, l1p0x, l1p0y, l1p1x, l1p1y );
}
#endif

static inline clpPoint64 clpGetActiveIntersectPoint( const clpActiveEdge * const restrict e1, const clpActiveEdge * const restrict e2 )
{
  double b1, b2;
  clpPoint64 pt;

  if( e1->slopedx == e2->slopedx )
    pt = e1->top;
  else if( e1->slopedx == 0.0 )
  {
    if( clpIsHorizontal( e2 ) )
      pt = clpBuildPoint64( e1->bot.x, e2->bot.y );
    else
    {
      b2 = e2->bot.y - ( e2->bot.x / e2->slopedx );
      pt = clpBuildPoint64( e1->bot.x, (int64_t)( round( ( e1->bot.x / e2->slopedx ) + b2 ) ) );
    }
  }
  else if( e2->slopedx == 0.0 )
  {
    if( clpIsHorizontal( e1 ) )
      pt = clpBuildPoint64( e2->bot.x, e1->bot.y );
    else
    {
      b1 = e1->bot.y - ( e1->bot.x / e1->slopedx );
      pt = clpBuildPoint64( e2->bot.x, (int64_t)( round( e2->bot.x / e1->slopedx + b1 ) ) );
    }
  }
  else
  {
#if 0
    double ptx, q;
    double e1dxabs, e2dxabs;
    e1dxabs = fabs( e1->slopedx );
    e2dxabs = fabs( e2->slopedx );
    b1 = e1->bot.x - ( e1->bot.y * e1->slopedx );
    b2 = e2->bot.x - ( e2->bot.y * e2->slopedx );
    q = ( b2 - b1 ) / ( e1->slopedx - e2->slopedx );
    /* Use the more vertical slope for improved accuracy */
    if( e1dxabs < e2dxabs )
      ptx = ( e1->slopedx * q ) + b1;
    else
      ptx = ( e2->slopedx * q ) + b2;
    pt = clpBuildPoint64( (int64_t)( round( ptx ) ), (int64_t)( round( q ) ) );
#else
    int64_t hitpt[2];
    clpVertexLineLineIntersection_C( hitpt, e1->bot.x, e1->bot.y, e1->top.x, e1->top.y, e2->bot.x, e2->bot.y, e2->top.x, e2->top.y );
    pt = clpBuildPoint64( hitpt[0], hitpt[1] );
#endif
  }

  return pt;
}


/*
static inline clpPoint64 clpGetClosestPointOnSegment( const clpPoint64 *offPt, const clpPoint64 *seg1, const clpPoint64 *seg2 )
{
  double q, sdx, sdy;
  if( ( seg1->x == seg2->x ) && ( seg1->y == seg2->y ) )
    return *seg1;
  sdx = (double)( seg2->x - seg1->x );
  sdy = (double)( seg2->y - seg1->y );
  q = ( (double)( offPt->x - seg1->x ) * (double)( seg2->x - seg1->x ) ) + ( (double)( offPt->y - seg1->y ) * (double)( seg2->y - seg1->y ) );
  q /= ( sdx * sdx ) + ( sdy * sdy );
  q = fmin( fmax( q, 0.0 ), 1.0 );
  sdx = ( ( 1.0 - q ) * (double)seg1->x ) + ( q * (double)seg2->x );
  sdy = ( ( 1.0 - q ) * (double)seg1->y ) + ( q * (double)seg2->y );
  return clpBuildPoint64( (int64_t)nearbyint( sdx ), (int64_t)nearbyint( sdy ) );
}
*/


static inline void clpAddNewIntersectNode( clpClipper *clp, clpActiveEdge *e1, clpActiveEdge *e2, int64_t scantopy )
{
  double e1dxabs, e2dxabs;
  clpPoint64 pt;
  clpIntersectNode *newnode;

  e1dxabs = fabs( e1->slopedx );
  e2dxabs = fabs( e2->slopedx );
  pt = clpGetActiveIntersectPoint( e1, e2 );

  /* Rounding errors can occasionally place the calculated intersection */
  /* point either below or above the scanbeam, so check and correct */
  if( pt.y > clp->scanbeambottomy )
  {
    /* Clamp to bottom of scanbeam */
    pt.y = clp->scanbeambottomy;
    /* Use the more vertical of the 2 edges to derive pt.x */
    if( e1dxabs < e2dxabs )
      pt.x = clpGetActiveTopX( e1, clp->scanbeambottomy );
    else
      pt.x = clpGetActiveTopX( e2, clp->scanbeambottomy );
  }
  else if( pt.y < scantopy )
  {
    /* Clamp to top of scanbeam */
    pt.y = scantopy;
    if( e1->top.y == scantopy )
      pt.x = e1->top.x;
    else if( e2->top.y == scantopy )
      pt.x = e2->top.x;
    else if( e1dxabs < e2dxabs )
      pt.x = e1->scanbeamx;
    else
      pt.x = e2->scanbeamx;
  }

  newnode = clpIntersectNodeVector_Access( &clp->intersectnodelist, clpIntersectNodeVector_Add( &clp->intersectnodelist ) );
  newnode->pt = pt;
  newnode->edge1 = e1;
  newnode->edge2 = e2;

  return;
}


/* TODO: This is the main loop where we could massively improve scalability by keeping an inactive queue of edges */
static inline int clpBuildIntersectList( clpClipper *clp, const int64_t scantopy )
{
  clpActiveEdge *left, *right, *leftend, *rightend, *curbase, *tmp;
  clpActiveEdge *prevbase;

  if( !clp->activeedgelist || !(clp->activeedgelist)->aelnext )
    return 0;

  /* Calculate edge positions at the top of the current scanbeam, and from this */
  /* we will determine the intersections required to reach these new positions */
  clpAdjustCurrXAndCopyToSEL( clp, scantopy );

  /* Find all edge intersections in the current scanbeam using a stable merge */
  /* sort that ensures only adjacent edges are intersecting. Intersect info is */
  /* stored in FIntersectList ready to be processed in ProcessIntersectList. */
  /* Re merge sorts see https://stackoverflow.com/a/46319131/359538 */

  for( left = clp->sortededgelist ; left && left->jump ; left = clp->sortededgelist )
  {
    for( prevbase = 0 ; left && left->jump ; prevbase = curbase, left = rightend )
    {
      curbase = left;
      right = left->jump;
      leftend = right;
      rightend = right->jump;
      left->jump = rightend;
      if( ( left == leftend ) || ( right == rightend ) )
        continue;
      for( ; ; )
      {
        /* Order is preserved, 'right' on right of 'left', nothing to see here */
        if( right->scanbeamx >= left->scanbeamx )
        {
          left = left->selnext;
          if( left == leftend )
            break;
        }
        else
        {
          /* If 'right' is on the left of 'left', we have intersections */
          /* Check all edges on the left of "right", until we hit "left" */
          for( tmp = right->selprev ; ; tmp = tmp->selprev )
          {
            /* Get intersection point between tmp and right */
            clpAddNewIntersectNode( clp, tmp, right, scantopy );
            if( tmp == left )
              break;
          }
          tmp = right;
          right = clpExtractFromSEL( tmp );
          leftend = right;
          clpInsert1Before2InSEL( tmp, left );
          if( left == curbase )
          {
            curbase = tmp;
            curbase->jump = rightend;
            if( !prevbase )
              clp->sortededgelist = curbase;
            else
              prevbase->jump = curbase;
          }
          /* 'right' was updated, so check for termination */
          if( right == rightend )
            break;
        }
      }
    }
  }
  return ( clpIntersectNodeVector_GetCount( &clp->intersectnodelist ) > 0 );
}

static inline void clpProcessIntersectList( clpClipper *clp )
{
  size_t nodeindex, nodecount;
  clpIntersectNode *node0, *node1, *nodeend;
  clpIntersectNode nodetmp;
  clpOutPt *op1, *op2;
  clpIntersectNode *nodelist;

  /* We now have a list of intersections required so that edges will be */
  /* correctly positioned at the top of the scanbeam. However, it's important */
  /* that edge intersections are processed from the bottom up, but it's also */
  /* crucial that intersections only occur between adjacent edges */

  nodecount = clpIntersectNodeVector_GetCount( &clp->intersectnodelist );
  nodelist = clpIntersectNodeVector_Access( &clp->intersectnodelist, 0 );

#if CLP_CONFIG_ENABLE_RADIX_SORT
  if( nodecount < CLP_CONFIG_RADIX_SORT_COUNT_THRESHOLD )
    clpIntersectNodeHybridSort( nodelist, 0, nodecount, 313893203 );
  else
  {
 #if CLP_CONFIG_SORT_INODE_X
    clpIntersectNodeRadixSortX( nodelist, 0, nodecount, 64 );
 #endif
    clpIntersectNodeRadixSortY( nodelist, 0, nodecount, 64 );
  }
#else
  clpIntersectNodeHybridSort( nodelist, 0, nodecount, 313893203 );
#endif

  /* Now as we process these intersections, we must sometimes adjust the order */
  /* to ensure that intersecting edges are always adjacent */

  nodecount = clpIntersectNodeVector_GetCount( &clp->intersectnodelist );
  nodeend = clpIntersectNodeVector_Access( &clp->intersectnodelist, nodecount );

  for( nodeindex = 0 ; nodeindex < nodecount ; nodeindex++ )
  {
    node0 = clpIntersectNodeVector_Access( &clp->intersectnodelist, nodeindex );
#if CLP_DEBUG_VERBOSE
    printf( "ProcessIntersectList %d %d\n", (int)node0->pt.x, (int)node0->pt.y );
#endif
    if( !clpEdgesAdjacentInAEL( node0 ) )
    {
      node1 = node0 + 1;
      while( ( node1 != nodeend ) && !clpEdgesAdjacentInAEL( node1 ) )
        ++node1;
      if( node1 != nodeend )
      {
        nodetmp = *node0;
        *node0 = *node1;
        *node1 = nodetmp;
      }
    }

    clpIntersectEdges( clp, node0->edge1, node0->edge2, node0->pt );
    clpSwapPositionsInAEL( clp, node0->edge1, node0->edge2 );

    if( clpTestJoinWithPrev2( node0->edge2, &node0->pt ) )
    {
      op1 = clpAddOutPt( clp, node0->edge2->aelprev, node0->pt );
      op2 = clpAddOutPt( clp, node0->edge2, node0->pt );
      if( op1 != op2 )
        clpAddJoin( clp, op1, op2 );
    }
    else if( clpTestJoinWithNext2( node0->edge1, &node0->pt ) )
    {
      op1 = clpAddOutPt( clp, node0->edge1, node0->pt );
      op2 = clpAddOutPt( clp, node0->edge1->aelnext, node0->pt );
      if( op1 != op2 )
        clpAddJoin( clp, op1, op2 );
    }
  }
  return;
}


static inline void clpDoIntersections( clpClipper *clp, const int64_t scantopy )
{
  if( clpBuildIntersectList( clp, scantopy ) )
  {
    clpProcessIntersectList( clp );
    clpIntersectNodeVector_Clear( &clp->intersectnodelist );
  }
  return;
}


static inline int clpHorzIsSpike( const clpActiveEdge *horzedge )
{
  clpPoint64 nextpt;
  nextpt = clpGetNextVertex( horzedge )->pt;
  return ( ( nextpt.y == horzedge->bot.y ) && ( ( horzedge->bot.x < horzedge->top.x ) != ( horzedge->top.x < nextpt.x ) ) );
}

static int clpTrimHorz( clpActiveEdge *horzedge, int preservecollinear )
{
  int result;
  clpPoint64 pt;
  result = 0;
  pt = clpGetNextVertex(horzedge)->pt;
  while( pt.y == horzedge->top.y )
  {
    /* Always trim 180 deg. spikes (in closed paths) */
    /* But otherwise break if preservecollinear = 1 */
    if( preservecollinear && ( ( pt.x < horzedge->top.x ) != ( horzedge->bot.x < horzedge->top.x ) ) )
      break;
    horzedge->vertextop = clpGetNextVertex( horzedge );
    horzedge->top = pt;
    result = 1;
    if( clpIsActiveEdgeMaxima( horzedge ) )
      break;
    pt = clpGetNextVertex(horzedge)->pt;
  }
  if( result )
    clpActiveUpdateDx( horzedge ); /* +/-infinity */
  return result;
}

static int clpResetHorzDirection( clpClipper *clp, const clpActiveEdge *horz, const clpActiveEdge *max_pair, int64_t *horz_left, int64_t *horz_right )
{
  clpActiveEdge *e;
  if( horz->bot.x == horz->top.x )
  {
    /* The horizontal edge is going nowhere */
    *horz_left = horz->scanbeamx;
    *horz_right = horz->scanbeamx;
    e = horz->aelnext;
    while( e && ( e != max_pair ) )
      e = e->aelnext;
    return ( e != 0 );
  }
  else if( horz->scanbeamx < horz->top.x )
  {
    *horz_left = horz->scanbeamx;
    *horz_right = horz->top.x;
    return 1;
  }
  else
  {
    *horz_left = horz->top.x;
    *horz_right = horz->scanbeamx;
    return 0;  /* Right to left */
  }
}


static void clpAddTrialHorzJoin( clpClipper *clp, clpOutPt *op )
{
  /* Make sure 'op' isn't added more than once */
  if( !( (op->outrec)->openflag ) && !clpOutPtInTrialHorzList( clp, op ) )
    clp->horzjoinerlist = clpNewJoiner( clp, op, 0, clp->horzjoinerlist );
  return;
}


  /*******************************************************************************
      * Notes: Horizontal edges (HEs) at scanline intersections (ie at the top or    *
      * bottom of a scanbeam) are processed as if layered.The order in which HEs     *
      * are processed doesn't matter. HEs intersect with the bottom vertices of      *
      * other HEs[#] and with non-horizontal edges [*]. Once these intersections     *
      * are completed, intermediate HEs are 'promoted' to the next edge in their     *
      * bounds, and they in turn may be intersected[%] by other HEs.                 *
      *                                                                              *
      * eg: 3 horizontals at a scanline:    /   |                     /           /  *
      *              |                     /    |     (HE3)o ========%========== o   *
      *              o ======= o(HE2)     /     |         /         /                *
      *          o ============#=========*======*========#=========o (HE1)           *
      *         /              |        /       |       /                            *
      *******************************************************************************/

static inline void clpDoHorizontal( clpClipper *clp, clpActiveEdge *horz )
{
  clpPoint64 pt;
  int horzisopen;
  int64_t y;
  clpVertex *vertex_max;
  clpActiveEdge *max_pair;
  int64_t horz_left, horz_right;
  int is_left_to_right;
  clpOutPt *op, *op2;
  clpActiveEdge *e;

  horzisopen = clpIsActiveEdgeOpen( horz );
  y = horz->bot.y;
  vertex_max = 0;
  max_pair = 0;

  if( !horzisopen )
  {
    vertex_max = clpGetCurrYMaximaVertex( horz );
    if( vertex_max )
    {
      max_pair = clpGetHorzMaximaPair( horz, vertex_max );
      /* Remove 180 deg.spikes and also simplify */
      /* Consecutive horizontals when preservecollinear = true */
      if( vertex_max != horz->vertextop )
        clpTrimHorz( horz, clp->preservecollinearflag );
    }
  }

  is_left_to_right = clpResetHorzDirection( clp, horz, max_pair, &horz_left, &horz_right );

  if( clpIsHotEdge( horz ) )
    clpAddOutPt( clp, horz, clpBuildPoint64( horz->scanbeamx, y ) );

  for( ; ; ) /* All consecutive horizontal edges */
  {
    if( horzisopen && clpIsActiveEdgeMaxima( horz ) && !clpIsActiveEdgeOpenEnd( horz ) )
    {
      vertex_max = clpGetCurrYMaximaVertex( horz );
      if( vertex_max )
        max_pair = clpGetHorzMaximaPair( horz, vertex_max );
    }

    if( is_left_to_right )
      e = horz->aelnext;
    else
      e = horz->aelprev;

    while( e )
    {
      if( e == max_pair )
      {
        if( clpIsHotEdge( horz ) )
        {
          while( horz->vertextop != e->vertextop )
          {
            clpAddOutPt( clp, horz, horz->top );
            clpUpdateEdgeIntoAEL( clp, horz );
          }
          if( is_left_to_right )
            op = clpAddLocalMaxPoly( clp, horz, e, horz->top );
          else
            op = clpAddLocalMaxPoly( clp, e, horz, horz->top );
          if( op && !clpIsActiveEdgeOpen( horz ) && clpPoint64Equal( &op->pt, &horz->top ) )
            clpAddTrialHorzJoin( clp, op );
        }
        clpDeleteFromAEL( clp, e );
        clpDeleteFromAEL( clp, horz );
        return;
      }

      /* If horzEdge is a maxima, keep going until we reach */
      /* its maxima pair, otherwise check for break conditions */
      if( ( vertex_max != horz->vertextop ) || clpIsActiveEdgeOpenEnd( horz ) )
      {
        /* Otherwise stop when activeedge is beyond the end of the horizontal line */
        if( ( is_left_to_right && ( e->scanbeamx > horz_right ) ) || ( !is_left_to_right && ( e->scanbeamx < horz_left ) ) )
          break;

        if( ( e->scanbeamx == horz->top.x ) && !clpIsHorizontal( e ) )
        {
          pt = clpGetNextVertex( horz )->pt;

          if( is_left_to_right )
          {
            /* With open paths we'll only break once past horz's end */
            if( clpIsActiveEdgeOpen( e ) && !clpIsSamePolyType( e, horz ) && !clpIsHotEdge( e ) )
            {
              if( clpGetActiveTopX( e, pt.y ) > pt.x )
                break;
            }
            /* Otherwise we'll only break when horz's outslope is greater than edge's */
            else if( clpGetActiveTopX( e, pt.y ) >= pt.x )
              break;
          }
          else
          {
            if( clpIsActiveEdgeOpen( e ) && !clpIsSamePolyType( e, horz ) && !clpIsHotEdge( e ) )
            {
              if( clpGetActiveTopX( e, pt.y ) < pt.x )
                break;
            }
            else if( clpGetActiveTopX( e, pt.y ) <= pt.x )
              break;
          }
        }
      }

      pt = clpBuildPoint64( e->scanbeamx, horz->bot.y );

      if( is_left_to_right )
      {
        if( clpIsActiveEdgeOpen( e ) && ( e->top.y == y ) )
          op = 0; /* Pass over the top of horz or maxpair open paths */
        else
          op = clpIntersectEdges( clp, horz, e, pt );
        clpSwapPositionsInAEL( clp, horz, e );

        if( clpIsHotEdge( horz ) && op && !clpIsActiveEdgeOpen( horz ) && clpPoint64Equal( &op->pt, &pt ) )
          clpAddTrialHorzJoin( clp, op );

        if( !clpIsHorizontal( e ) && clpTestJoinWithPrev1( e, y ) )
        {
          op = clpAddOutPt( clp, e->aelprev, pt );
          op2 = clpAddOutPt( clp, e, pt );
          clpAddJoin( clp, op, op2 );
        }

        horz->scanbeamx = e->scanbeamx;
        e = horz->aelnext;
      }
      else
      {
        if( clpIsActiveEdgeOpen( e ) && ( e->top.y == y ) )
          op = 0; /* Pass over the top of horz or maxpair open paths */
        else
          op = clpIntersectEdges( clp, e, horz, pt );
        clpSwapPositionsInAEL( clp, e, horz );

        if( clpIsHotEdge( horz ) && op && !clpIsActiveEdgeOpen( horz ) && clpPoint64Equal( &op->pt, &pt ) )
          clpAddTrialHorzJoin( clp, op );

        if( !clpIsHorizontal( e ) && clpTestJoinWithNext1( e, y ) )
        {
          op = clpAddOutPt( clp, e, pt );
          op2 = clpAddOutPt( clp, e->aelnext, pt );
          clpAddJoin( clp, op, op2 );
        }

        horz->scanbeamx = e->scanbeamx;
        e = horz->aelprev;
      }
    }

    /* Check if we've finished with (consecutive) horizontals */
    if( horzisopen && clpIsActiveEdgeOpenEnd( horz ) ) /* Open at top */
    {
      if( clpIsHotEdge( horz ) )
      {
        clpAddOutPt( clp, horz, horz->top );
        if( clpIsFront( horz ) )
          horz->outrec->front_edge = 0;
        else
          horz->outrec->back_edge = 0;
        horz->outrec = 0;
      }
      clpDeleteFromAEL( clp, horz );
      return;
    }
    else if( clpGetNextVertex( horz )->pt.y != horz->top.y )
      break;

    /* Still more horizontals in bound to process */
    if( clpIsHotEdge( horz ) )
      clpAddOutPt( clp, horz, horz->top );
    clpUpdateEdgeIntoAEL( clp, horz );

    if( clp->preservecollinearflag && clpHorzIsSpike( horz ) )
      clpTrimHorz( horz, 1 );

    is_left_to_right = clpResetHorzDirection( clp, horz, max_pair, &horz_left, &horz_right );
  }

  if( clpIsHotEdge( horz ) )
  {
    op = clpAddOutPt( clp, horz, horz->top );
    if( !clpIsActiveEdgeOpen( horz ) )
      clpAddTrialHorzJoin( clp, op );
  }
  else
    op = 0;

  if( ( horzisopen && !clpIsActiveEdgeOpenEnd( horz ) ) || ( !horzisopen && ( vertex_max != horz->vertextop ) ) )
  {
    clpUpdateEdgeIntoAEL( clp, horz ); /* This is the end of an intermediate horz */
    if( clpIsActiveEdgeOpen( horz ) )
      return;

    if( is_left_to_right && clpTestJoinWithNext1( horz, y ) )
    {
      op2 = clpAddOutPt( clp, horz->aelnext, horz->bot );
      clpAddJoin( clp, op, op2 );
    }
    else if( !is_left_to_right && clpTestJoinWithPrev1( horz, y ) )
    {
      op2 = clpAddOutPt( clp, horz->aelprev, horz->bot );
      clpAddJoin( clp, op2, op );
    }
  }
  else if( clpIsHotEdge( horz ) )
    clpAddLocalMaxPoly( clp, horz, max_pair, horz->top );
  else
  {
    clpDeleteFromAEL( clp, max_pair );
    clpDeleteFromAEL( clp, horz );
  }
  return;
}

static void clpDoHorizontalPopAll( clpClipper *clp )
{
  clpActiveEdge *ae;
  while( ( ae = clpPopHorz( clp ) ) )
  {
#if CLP_DEBUG_VERBOSE
    printf( "  DoHorizontal %d,%d %d,%d ~ %d %f %d %d %d\n", (int)ae->bot.x, (int)ae->bot.y, (int)ae->top.x, (int)ae->top.y, (int)ae->scanbeamx, ae->slopedx, ae->windingdelta, ae->windingcount, ae->windingcount2 );
#endif
    clpDoHorizontal( clp, ae );
  }
  return;
}


static inline clpActiveEdge *clpDoMaxima( clpClipper *clp, clpActiveEdge *e )
{
  clpActiveEdge *next_e, *prev_e, *max_pair;
  prev_e = e->aelprev;
  next_e = e->aelnext;
  if( clpIsActiveEdgeOpenEnd( e ) )
  {
    if( clpIsHotEdge( e ) )
      clpAddOutPt( clp, e, e->top );
    if( !clpIsHorizontal( e ) )
    {
      if( clpIsHotEdge( e ) )
      {
        if( clpIsFront( e ) )
          e->outrec->front_edge = 0;
        else
          e->outrec->back_edge = 0;
        e->outrec = 0;
      }
      clpDeleteFromAEL( clp, e );
    }
    return next_e;
  }
  else
  {
    max_pair = clpGetMaximaPair( e );
    if( !max_pair )
      return next_e;
  }

  /* Only non-horizontal maxima here */
  /* Process any edges between maxima pair */
  while( next_e != max_pair )
  {
    clpIntersectEdges( clp, e, next_e, e->top );
    clpSwapPositionsInAEL( clp, e, next_e );
    next_e = e->aelnext;
  }

  if( clpIsActiveEdgeOpen( e ) )
  {
    if( clpIsHotEdge( e ) )
      clpAddLocalMaxPoly( clp, e, max_pair, e->top );
    clpDeleteFromAEL( clp, max_pair );
    clpDeleteFromAEL( clp, e );
    return ( prev_e ? prev_e->aelnext : clp->activeedgelist );
  }

  if( clpIsHotEdge( e ) )
    clpAddLocalMaxPoly( clp, e, max_pair, e->top );

  clpDeleteFromAEL( clp, e );
  clpDeleteFromAEL( clp, max_pair );
  return ( prev_e ? prev_e->aelnext : clp->activeedgelist );
}


static inline void clpDoTopOfScanbeam( clpClipper *clp, const int64_t y )
{
  clpActiveEdge *e;
  /* sortededgelist is reused to flag horizontals (see clpPushHorz below) */
  clp->sortededgelist = 0;
  e = clp->activeedgelist;
  while( e )
  {
    /* Active edge will never be horizontal here */
    if( e->top.y == y )
    {
      e->scanbeamx = e->top.x;
      if( clpIsActiveEdgeMaxima( e ) )
      {
        e = clpDoMaxima( clp, e );
        continue;
      }
      else
      {
        /* Intermediate vertex */
        if( clpIsHotEdge( e ) )
          clpAddOutPt( clp, e, e->top );
        clpUpdateEdgeIntoAEL( clp, e );
        if( clpIsHorizontal( e ) )
          clpPushHorz( clp, e );  /* Horizontals are processed later */
      }
    }
    else /* It wasn't the top of the edge */
      e->scanbeamx = clpGetActiveTopX( e, y );
    e = e->aelnext;
  }
  return;
}


static int clpGetHorzExtendedHorzSeg( clpOutPt **inoutop, clpOutPt **inoutop2 )
{
  int retval;
  clpOutPt *op, *op2;
  clpOutRec *outrec;
  op = *inoutop;
  outrec = clpGetRealOutRec( op->outrec );
  op2 = op;
  if( outrec->front_edge )
  {
    while( ( op->prev != outrec->pts ) && ( (op->prev)->pt.y == op->pt.y ) )
      op = op->prev;
    while( ( op2 != outrec->pts ) && ( (op2->next)->pt.y == op2->pt.y ) )
      op2 = op2->next;
    retval = ( op2 != op );
  }
  else
  {
    while( ( op->prev != op2 ) && ( (op->prev)->pt.y == op->pt.y ) )
      op = op->prev;
    while( ( op2->next != op ) && ( (op2->next)->pt.y == op2->pt.y ) )
      op2 = op2->next;
    retval = ( ( op2 != op ) && ( op2->next != op ) );
  }
  *inoutop = op;
  *inoutop2 = op2;
  return retval;
}


static inline int clpHorzEdgesOverlap( int64_t x1a, int64_t x1b, int64_t x2a, int64_t x2b )
{
  const int64_t minOverlap = 2;
  if( x1a > ( x1b + minOverlap ) )
  {
    if( x2a > ( x2b + minOverlap ) )
      return !( ( x1a <= x2b ) || ( x2a <= x1b ) );
    else
      return !( ( x1a <= x2a ) || ( x2b <= x1b ) );
  }
  else if( x1b > ( x1a + minOverlap ) )
  {
    if( x2a > ( x2b + minOverlap ) )
      return !( ( x1b <= x2b ) || ( x2a <= x1a ) );
    else
      return !( ( x1b <= x2a ) || ( x2b <= x1a ) );
  }
  else
    return 0;
}


static inline int clpValueBetween( int64_t val, int64_t end1, int64_t end2 )
{
  /* Accommodates axis aligned between where end1 == end2 */
  return ( ( ( val != end1 ) == ( val != end2 ) ) && ( ( val > end1 ) == ( val < end2 ) ) );
}

static inline int clpValueEqualOrBetween( int64_t val, int64_t end1, int64_t end2 )
{
  return ( ( val == end1 ) || ( val == end2 ) || ( ( val > end1 ) == ( val < end2 ) ) );
}

static inline int clpPointBetween( clpPoint64 pt, clpPoint64 corner1, clpPoint64 corner2 )
{
  /* Points may not be collinear */
  return clpValueEqualOrBetween( pt.x, corner1.x, corner2.x ) && clpValueEqualOrBetween( pt.y, corner1.y, corner2.y );
}


static inline void clpConvertHorzTrialsToJoins( clpClipper *clp )
{
  clpOutPt *op1a;
  clpOutPt *op1b;
  clpOutPt *op2a, *op2b;
  clpJoiner *joiner, *joinerParent;
  int joined;

  while( clp->horzjoinerlist )
  {
    joiner = clp->horzjoinerlist;
    clp->horzjoinerlist = joiner->nextH;
    op1a = joiner->op1;
    if( op1a->joiner == joiner )
      op1a->joiner = joiner->next1;
    else
    {
      joinerParent = clpFindJoinParent( joiner, op1a );
      if( joinerParent->op1 == op1a )
        joinerParent->next1 = joiner->next1;
      else
        joinerParent->next2 = joiner->next1;
    }
    clpFreeJoiner( clp, joiner );

    if( !clpGetHorzExtendedHorzSeg( &op1a, &op1b ) )
    {
      clpCleanCollinear( clp, op1a->outrec );
      continue;
    }

#if CLP_DEBUG_VERBOSE
    DebugDebugPrintPtChain( op1a, "ConvertHorzTrialsToJoins" );
#endif

    joined = 0;
    joiner = clp->horzjoinerlist;
    while( joiner )
    {
      op2a = joiner->op1;
      if( clpGetHorzExtendedHorzSeg( &op2a, &op2b ) && clpHorzEdgesOverlap( op1a->pt.x, op1b->pt.x, op2a->pt.x, op2b->pt.x ) )
      {
        /* Overlap found so promote to a 'real' join */
        joined = 1;
        if( clpPoint64Equal( &op1a->pt, &op2b->pt ) )
          clpAddJoin( clp, op1a, op2b );
        else if( clpPoint64Equal( &op1b->pt, &op2a->pt ) )
          clpAddJoin( clp, op1b, op2a );
        else if( clpPoint64Equal( &op1a->pt, &op2a->pt ) )
          clpAddJoin( clp, op1a, op2a );
        else if( clpPoint64Equal( &op1b->pt, &op2b->pt ) )
          clpAddJoin( clp, op1b, op2b );
        else if( clpValueBetween( op1a->pt.x, op2a->pt.x, op2b->pt.x ) )
          clpAddJoin( clp, op1a, clpInsertOp( clp, &op1a->pt, op2a ) );
        else if( clpValueBetween( op1b->pt.x, op2a->pt.x, op2b->pt.x ) )
          clpAddJoin( clp, op1b, clpInsertOp( clp, &op1b->pt, op2a ) );
        else if( clpValueBetween( op2a->pt.x, op1a->pt.x, op1b->pt.x ) )
          clpAddJoin( clp, op2a, clpInsertOp( clp, &op2a->pt, op1a ) );
        else if( clpValueBetween( op2b->pt.x, op1a->pt.x, op1b->pt.x ) )
          clpAddJoin( clp, op2b, clpInsertOp( clp, &op2b->pt, op1a ) );
        break;
      }
      joiner = joiner->nextH;
    }

#if CLP_DEBUG_VERBOSE
    if( !joined )
      printf( "ConvertHorzTrialsToJoins, Enter CleanCollinear()\n" );
#endif

    if( !joined )
      clpCleanCollinear( clp, op1a->outrec );
  }
  return;
}


static int clpCheckDisposeAdjacent( clpClipper *clp, clpOutPt **inoutop, const clpOutPt *guard, clpOutRec *outRec )
{
  clpOutPt *op;
  int result;
  op = *inoutop;
  result = 0;
  while( op->prev != op )
  {
#if CLP_DEBUG_VERBOSE
    printf( "  CheckDisposeAdjacent A %d %d\n", (int)op->pt.x, (int)op->pt.y );
#endif
    if( clpPoint64Equal( &op->pt, &(op->prev)->pt ) && ( op != guard ) && op->prev->joiner && !op->joiner )
    {
#if CLP_DEBUG_VERBOSE
      printf( "  CheckDisposeAdjacent AAA\n", (int)op->pt.x, (int)op->pt.y );
#endif
      if( op == outRec->pts )
        outRec->pts = op->prev;
      op = clpDisposeOutPt( clp, op );
      op = op->prev;
    }
    else
      break;
  }

  while( op->next != op )
  {
#if CLP_DEBUG_VERBOSE
    printf( "  CheckDisposeAdjacent B %d %d\n", (int)op->pt.x, (int)op->pt.y );
#endif
    if( clpPoint64Equal( &op->pt, &(op->next)->pt ) && ( op != guard ) && op->next->joiner && !op->joiner )
    {
#if CLP_DEBUG_VERBOSE
      printf( "  CheckDisposeAdjacent CCC\n", (int)op->pt.x, (int)op->pt.y );
#endif
      if( op == outRec->pts )
        outRec->pts = op->prev;
      op = clpDisposeOutPt( clp, op );
      op = op->prev;
    }
    else
      break;
  }
  *inoutop = op;
  return result;
}


static inline int clpIsValidPath( clpOutPt *op )
{
  return ( op && ( op->next != op ) );
}

static int clpCollinearSegsOverlap( const clpPoint64 *seg1a, const clpPoint64 *seg1b, const clpPoint64 *seg2a, const clpPoint64 *seg2b )
{
  /* Precondition: seg1 and seg2 are collinear */
  if( seg1a->x == seg1b->x )
  {
    if( ( seg2a->x != seg1a->x ) || ( seg2a->x != seg2b->x ) )
      return 0;
  }
  else if( seg1a->x < seg1b->x )
  {
    if( seg2a->x < seg2b->x )
    {
      if( ( seg2a->x >= seg1b->x ) || ( seg2b->x <= seg1a->x ) )
        return 0;
    }
    else
    {
      if( ( seg2b->x >= seg1b->x ) || ( seg2a->x <= seg1a->x ) )
        return 0;
    }
  }
  else
  {
    if( seg2a->x < seg2b->x )
    {
      if( ( seg2a->x >= seg1a->x ) || ( seg2b->x <= seg1b->x ) )
        return 0;
    }
    else
    {
      if( ( seg2b->x >= seg1a->x ) || ( seg2a->x <= seg1b->x ) )
        return 0;
    }
  }

  if( seg1a->y == seg1b->y )
  {
    if( ( seg2a->y != seg1a->y ) || ( seg2a->y != seg2b->y ) )
      return 0;
  }
  else if( seg1a->y < seg1b->y )
  {
    if( seg2a->y < seg2b->y )
    {
      if( ( seg2a->y >= seg1b->y ) || ( seg2b->y <= seg1a->y ) )
        return 0;
    }
    else
    {
      if( ( seg2b->y >= seg1b->y ) || ( seg2a->y <= seg1a->y ) )
        return 0;
    }
  }
  else
  {
    if( seg2a->y < seg2b->y )
    {
      if( ( seg2a->y >= seg1a->y ) || ( seg2b->y <= seg1b->y ) )
        return 0;
    }
    else
    {
      if( ( seg2b->y >= seg1a->y ) || ( seg2a->y <= seg1b->y ) )
        return 0;
    }
  }
  return 1;
}


/* Perpendicular distance of point to line */
static inline double clpDistanceFromLineSqrd( const clpPoint64 *pt, const clpPoint64 *ln1, const clpPoint64 *ln2 )
{
  double A, B, C;
  A = (double)( ln1->y - ln2->y );
  B = (double)( ln2->x - ln1->x );
  C = ( A * ln1->x ) + ( B * ln1->y );
  C = ( A * pt->x ) + ( B * pt->y ) - C;
  return ( C * C ) / ( ( A * A ) + ( B * B ) );
}


static inline clpOutRec *clpProcessJoin( clpClipper *clp, clpJoiner *joiner )
{
  clpOutPt *op1, *op2;
  clpOutRec *or1;
  clpOutRec *or2;
  clpOutPt *opA, *opB;
  clpOutRec *result;

  op1 = joiner->op1;
  op2 = joiner->op2;
  or1 = clpGetRealOutRec( op1->outrec );
  or2 = clpGetRealOutRec( op2->outrec );
  clpDeleteJoin( clp, joiner );

#if CLP_DEBUG_VERBOSE
  DebugDebugPrintPtChain( op1, "  ProcessJoin Op1" );
  DebugDebugPrintPtChain( op2, "  ProcessJoin Op2" );
#endif

#if CLP_DEBUG_VERBOSE
  printf( "  ProcessJoin\n" );
#endif

  if( or2->pts == 0 )
  {
#if CLP_DEBUG_VERBOSE
    printf( "  ProcessJoin A\n" );
#endif
    return or1;
  }
  else if( !clpIsValidClosedPath( op2 ) )
  {
#if CLP_DEBUG_VERBOSE
    printf( "  ProcessJoin B\n" );
#endif
    clpCleanCollinear( clp, or2 );
    return or1;
  }
  else if( ( or1->pts == 0 ) || !clpIsValidClosedPath( op1 ) )
  {
#if CLP_DEBUG_VERBOSE
    printf( "  ProcessJoin C\n" );
#endif
    clpCleanCollinear( clp, or1 );
    return or2;
  }
  else if( ( or1 == or2 ) && ( ( op1 == op2 ) || ( op1->next == op2 ) || ( op1->prev == op2 ) ) )
  {
#if CLP_DEBUG_VERBOSE
    printf( "  ProcessJoin D\n" );
#endif
    return or1;
  }

#if CLP_DEBUG_VERBOSE
  DebugDebugPrintPtChain( op1, "  ProcessJoin Op1 A" );
  DebugDebugPrintPtChain( op2, "  ProcessJoin Op2 A" );
#endif

  clpCheckDisposeAdjacent( clp, &op1, op2, or1 );
  clpCheckDisposeAdjacent( clp, &op2, op1, or2 );

#if CLP_DEBUG_VERBOSE
  DebugDebugPrintPtChain( op1, "  ProcessJoin Op1 B" );
  DebugDebugPrintPtChain( op2, "  ProcessJoin Op2 B" );
#endif

  if( ( op1->next == op2 ) || ( op2->next == op1 ) )
  {
#if CLP_DEBUG_VERBOSE
    printf( "  ProcessJoin E\n" );
#endif
    return or1;
  }

  result = or1;
  for( ; ; )
  {
#if CLP_DEBUG_VERBOSE
    DebugDebugPrintPtChain( result->pts, "  ProcessJoin TMP000" );
#endif

#if CLP_DEBUG_VERBOSE
    printf( "  ProcessJoin %d,%d %d,%d\n", (int)op1->pt.x, (int)op1->pt.y, (int)op2->pt.x, (int)op2->pt.y );
#endif

    if( !clpIsValidPath( op1 ) || !clpIsValidPath( op2 ) || ( ( or1 == or2 ) && ( ( op1->prev == op2 ) || ( op1->next == op2 ) ) ) )
      return or1;
    if( clpPoint64Equal( &(op1->prev)->pt, &(op2->next)->pt ) || ( ( clpCrossProduct( &(op1->prev)->pt, &op1->pt, &(op2->next)->pt ) == 0 ) && clpCollinearSegsOverlap( &(op1->prev)->pt, &op1->pt, &op2->pt, &(op2->next)->pt ) ) )
    {
      if( or1 == or2 )
      {
        /* Split required */
        /* make sure op1.prev and op2.next match positions */
        /* by inserting an extra vertex if needed */
        if( !clpPoint64Equal( &(op1->prev)->pt, &(op2->next)->pt ) )
        {
#if CLP_DEBUG_VERBOSE
          printf( "  ProcessJoin AAA\n" );
#endif
          if( clpPointBetween( (op1->prev)->pt, op2->pt, (op2->next)->pt ) )
            op2->next = clpInsertOp( clp, &(op1->prev)->pt, op2 );
          else
            op1->prev = clpInsertOp( clp, &(op2->next)->pt, op1->prev );
        }

        /* current              to     new           */
        /* op1.p[opA] >>> op1   ...    opA \   / op1 */
        /* op2.n[opB] <<< op2   ...    opB /   \ op2 */

        opA = op1->prev;
        opB = op2->next;
        opA->next = opB;
        opB->prev = opA;
        op1->prev = op2;
        op2->next = op1;
        clpCompleteSplit( clp, op1, opA, or1 );
#if CLP_DEBUG_VERBOSE
        printf( "  ProcessJoin BBB\n" );
#endif
      }
      else
      {
        /* Join, not split */
        opA = op1->prev;
        opB = op2->next;
        opA->next = opB;
        opB->prev = opA;
        op1->prev = op2;
        op2->next = op1;
        if( or1->idx < or2->idx )
        {
          or1->pts = op1;
          or2->pts = 0;
          if( or1->owner && ( !or2->owner || ( (or2->owner)->idx < (or1->owner)->idx ) ) )
            or1->owner = or2->owner;
          or2->owner = or1;
#if CLP_DEBUG_VERBOSE
          printf( "  ProcessJoin CCC\n" );
#endif
        }
        else
        {
          result = or2;
          or2->pts = op1;
          or1->pts = 0;
          if( or2->owner && ( !or1->owner || ( (or1->owner)->idx < (or2->owner)->idx ) ) )
            or2->owner = or1->owner;
          or1->owner = or2;
#if CLP_DEBUG_VERBOSE
          printf( "  ProcessJoin DDD\n" );
#endif
        }
      }
      break;
    }
    else if( clpPoint64Equal( &(op1->next)->pt, &(op2->prev)->pt ) || ( ( clpCrossProduct( &(op1->next)->pt, &op2->pt, &(op2->prev)->pt ) == 0 ) && clpCollinearSegsOverlap( &(op1->next)->pt, &op1->pt, &op2->pt, &(op2->prev)->pt ) ) )
    {
#if CLP_DEBUG_VERBOSE
      printf( "  ProcessJoin EEE\n" );
#endif
      if( or1 == or2 )
      {
        /* Split required */
        /* Make sure op2.prev and op1.next match positions */
        /* by inserting an extra vertex if needed */
#if CLP_DEBUG_VERBOSE
        printf( "  ProcessJoin FFF\n" );
#endif
        if( !clpPoint64Equal( &(op2->prev)->pt, &(op1->next)->pt ) )
        {
#if CLP_DEBUG_VERBOSE
          printf( "  ProcessJoin GGG\n" );
#endif
          if( clpPointBetween( (op2->prev)->pt, op1->pt, (op1->next)->pt ) )
            op1->next = clpInsertOp( clp, &(op2->prev)->pt, op1 );
          else
            op2->prev = clpInsertOp( clp, &(op1->next)->pt, op2->prev );
        }

        /* current              to     new           */
        /* op2.p[opA] >>> op2   ...    opA \   / op2 */
        /* op1.n[opB] <<< op1   ...    opB /   \ op1 */
        opA = op2->prev;
        opB = op1->next;
        opA->next = opB;
        opB->prev = opA;
        op2->prev = op1;
        op1->next = op2;
        clpCompleteSplit( clp, op1, opA, or1 );
      }
      else
      {
#if CLP_DEBUG_VERBOSE
        printf( "  ProcessJoin HHH\n" );
#endif
        /* Join, not split */
        opA = op1->next;
        opB = op2->prev;
        opA->prev = opB;
        opB->next = opA;
        op1->next = op2;
        op2->prev = op1;
        if( or1->idx < or2->idx )
        {
#if CLP_DEBUG_VERBOSE
          printf( "  ProcessJoin III\n" );
#endif
          or1->pts = op1;
          or2->pts = 0;
          if( or1->owner && ( !or2->owner || ( (or2->owner)->idx < (or1->owner)->idx ) ) )
            or1->owner = or2->owner;
          or2->owner = or1;
        }
        else
        {
#if CLP_DEBUG_VERBOSE
          printf( "  ProcessJoin JJJ\n" );
#endif
          result = or2;
          or2->pts = op1;
          or1->pts = 0;
          if( or2->owner && ( !or1->owner || ( (or1->owner)->idx < (or2->owner)->idx ) ) )
            or2->owner = or1->owner; 
          or1->owner = or2;
        }
      }
      break;
    }
    else if( clpPointBetween( (op1->next)->pt, op2->pt, (op2->prev)->pt ) && clpDistanceFromLineSqrd( &(op1->next)->pt, &op2->pt, &(op2->prev)->pt ) < 2.01 )
    {
#if CLP_DEBUG_VERBOSE
      printf( "  ProcessJoin KKK\n" );
#endif
      clpInsertOp( clp, &(op1->next)->pt, op2->prev );
      continue;
    }
    else if( clpPointBetween( (op2->next)->pt, op1->pt, (op1->prev)->pt ) && clpDistanceFromLineSqrd( &(op2->next)->pt, &op1->pt, &(op1->prev)->pt ) < 2.01 )
    {
#if CLP_DEBUG_VERBOSE
      printf( "  ProcessJoin LLL\n" );
#endif
      clpInsertOp( clp, &(op2->next)->pt, op1->prev );
      continue;
    }
    else if( clpPointBetween( (op1->prev)->pt, op2->pt, (op2->next)->pt ) && clpDistanceFromLineSqrd( &(op1->prev)->pt, &op2->pt, &(op2->next)->pt ) < 2.01 )
    {
#if CLP_DEBUG_VERBOSE
      printf( "  ProcessJoin MMM\n" );
#endif
      clpInsertOp( clp, &(op1->prev)->pt, op2 );
      continue;
    }
    else if( clpPointBetween( (op2->prev)->pt, op1->pt, (op1->next)->pt ) && clpDistanceFromLineSqrd( &(op2->prev)->pt, &op1->pt, &(op1->next)->pt ) < 2.01 )
    {
#if CLP_DEBUG_VERBOSE
      printf( "  ProcessJoin NNN\n" );
#endif
      clpInsertOp( clp, &(op2->prev)->pt, op1 );
      continue;
    }

    /* Something odd needs tidying up */
    if( clpCheckDisposeAdjacent( clp, &op1, op2, or1 ) )
      continue;
    else if( clpCheckDisposeAdjacent( clp, &op2, op1, or1 ) )
      continue;
    else if( !clpPoint64Equal( &(op1->prev)->pt, &(op2->next)->pt ) && ( clpDistanceSqr( &(op1->prev)->pt, &(op2->next)->pt ) < 2.01 ) )
    {
#if CLP_DEBUG_VERBOSE
      printf( "  ProcessJoin OOO\n" );
#endif
      op1->prev->pt = op2->next->pt;
      continue;
    }
    else if( !clpPoint64Equal( &(op1->next)->pt, &(op2->prev)->pt ) && ( clpDistanceSqr( &(op1->next)->pt, &(op2->prev)->pt ) < 2.01 ) )
    {
#if CLP_DEBUG_VERBOSE
      printf( "  ProcessJoin PPP\n" );
#endif
      op2->prev->pt = op1->next->pt;
      continue;
    }
    else
    {
#if CLP_DEBUG_VERBOSE
      printf( "  ProcessJoin QQQ\n" );
#endif
      /* OK, there doesn't seem to be a way to join after all */
      /* so just tidy up the polygons */
      or1->pts = op1;
      if( or2 != or1 )
      {
        or2->pts = op2;
        clpCleanCollinear( clp, or2 );
      }
      break;
    }
  }

#if CLP_DEBUG_VERBOSE
  DebugDebugPrintPtChain( result->pts, "  ProcessJoin Result" );
#endif

  return result;
}


static inline void clpProcessJoinerList( clpClipper *clp )
{
  size_t joinerindex, joinercount;
  clpOutRec *outrec;
  clpJoiner *j;
  joinercount = clpJoinerVector_GetCount( &clp->joinerlist );
  for( joinerindex = 0 ; joinerindex < joinercount ; joinerindex++ )
  {
    j = clpJoinerVector_Read( &clp->joinerlist, joinerindex );
    if( !j )
      continue;
    else if( !clp->succeessflag )
      clpFreeJoiner( clp, j );
    else
    {
      outrec = clpProcessJoin( clp, j );
      clpCleanCollinear( clp, outrec );
    }
  }
  clpJoinerVector_Clear( &clp->joinerlist );
  return;
}


static inline int clpExecuteInternal( clpClipper *clp, int cliptype, int fillrule, int use_polytrees )
{
  int64_t y;

  clp->cliptype = cliptype;
  clp->fillrule = fillrule;
  clp->polytreeflag = use_polytrees;
  clpReset( clp );

  if( cliptype == CLP_CLIPTYPE_NONE )
    return 1;
  if( !clpPopScanline( clp, &y ) )
    return 1;

  while( clp->succeessflag )
  {
#if CLP_DEBUG_VERBOSE
    printf( "ExecuteInternal, InsertLocalMinimaIntoAEL %d\n", (int)y );
#endif
    clpInsertLocalMinimaIntoAEL( clp, y );
    clpDoHorizontalPopAll( clp );
    if( clp->horzjoinerlist )
      clpConvertHorzTrialsToJoins( clp );
    clp->scanbeambottomy = y;  /* scanbeambottomy == bottom of scanbeam */
    if( !clpPopScanline( clp, &y ) )
      break;  /* y new top of scanbeam */
    clpDoIntersections( clp, y );
    clpDoTopOfScanbeam( clp, y );
    clpDoHorizontalPopAll( clp );
  }
  clpProcessJoinerList( clp );
  return clp->succeessflag;
}


////


static inline int clpBuildPath( clpOutPt *op, int reverseflag, int openflag, clpPath *path )
{
  size_t vertexcount, vertexalloc;
  clpPoint64 lastPt;
  clpOutPt *op2;

  path->vertexcount = 0;
  path->vertexlist = 0;
  if( ( op->next == op ) || ( !openflag && ( op->next == op->prev ) ) )
    return 0;
  vertexcount = 0;
  vertexalloc = 256;
  path->vertexlist = (int64_t *)malloc( ( vertexalloc + CLP_CONFIG_PATHVLIST_EXTRA_ALLOC ) * 2 * sizeof(int64_t) );
  if( reverseflag )
  {
    lastPt = op->pt;
    path->vertexlist[0] = lastPt.x;
    path->vertexlist[1] = lastPt.y;
    vertexcount++;
    for( op2 = op->prev ; op2 != op ; op2 = op2->prev )
    {
#if CLP_DEBUG_VERBOSE
      printf( "BuildPath %d %d\n", (int)op2->pt.x, (int)op2->pt.y );
#endif
      if( !clpPoint64Equal( &op2->pt, &lastPt ) )
      {
        lastPt = op2->pt;
        if( vertexcount >= vertexalloc )
        {
          /* Grow alloc by 25% plus 256 */
          vertexalloc = ( ( vertexalloc * 5 ) >> 2 ) + 256;
          path->vertexlist = (int64_t *)realloc( path->vertexlist, ( vertexalloc + CLP_CONFIG_PATHVLIST_EXTRA_ALLOC ) * 2 * sizeof(int64_t) );
        }
        path->vertexlist[(2*vertexcount)+0] = lastPt.x;
        path->vertexlist[(2*vertexcount)+1] = lastPt.y;
        vertexcount++;
      }
    }
  }
  else
  {
    op = op->next;
    lastPt = op->pt;
    path->vertexlist[0] = lastPt.x;
    path->vertexlist[1] = lastPt.y;
    vertexcount++;
    for( op2 = op->next ; op2 != op ; op2 = op2->next )
    {
#if CLP_DEBUG_VERBOSE
      printf( "BuildPath %d %d\n", (int)op2->pt.x, (int)op2->pt.y );
#endif
      if( !clpPoint64Equal( &op2->pt, &lastPt ) )
      {
        lastPt = op2->pt;
        if( vertexcount >= vertexalloc )
        {
          /* Grow alloc by 25% plus 256 */
          vertexalloc = ( ( vertexalloc * 5 ) >> 2 ) + 256;
          path->vertexlist = (int64_t *)realloc( path->vertexlist, ( vertexalloc + CLP_CONFIG_PATHVLIST_EXTRA_ALLOC ) * 2 * sizeof(int64_t) );
        }
        path->vertexlist[(2*vertexcount)+0] = lastPt.x;
        path->vertexlist[(2*vertexcount)+1] = lastPt.y;
        vertexcount++;
      }
    }
  }
  path->vertexcount = vertexcount;
  return 1;
}


static inline void clpBuildShape( clpClipper *clp, clpShape *shapeclosed, clpShape *shapeopen )
{
  size_t closedindex, openindex;
  size_t outrecindex, outreclistcount;
  clpOutRec *outrec;

  outreclistcount = clpOutRecVector_GetCount( &clp->outreclist );
  closedindex = 0;
  openindex = 0;
  shapeclosed->pathcount = 0;
  shapeclosed->pathlist = (clpPath *)malloc( outreclistcount * sizeof(clpPath) );
  if( shapeopen )
  {
    shapeopen->pathcount = 0;
    shapeopen->pathlist = (clpPath *)malloc( outreclistcount * sizeof(clpPath) );
  }

  for( outrecindex = 0 ; outrecindex < outreclistcount ; outrecindex++ )
  {
    outrec = clpOutRecVector_Read( &clp->outreclist, outrecindex );
    if( outrec->pts == 0 )
      continue;
    if( shapeopen && outrec->openflag )
    {
      if( clpBuildPath( outrec->pts, clp->reversesolutionflag, 1, &shapeopen->pathlist[openindex] ) )
        openindex++;
    }
    else
    {
      /* Closed paths should always return a positive orientation */
      if( clpBuildPath( outrec->pts, ( clp->reversesolutionflag != clp->orientationreversedflag ), 0, &shapeclosed->pathlist[closedindex] ) )
        closedindex++;
    }
  }
  shapeclosed->pathcount = closedindex;
  if( shapeopen )
    shapeopen->pathcount = openindex;
  return;
}


////


#if CLP_CONFIG_PATHVLIST_EXTRA_ALLOC < 4
 #error pwPointInContour64_Padding() requires CLP_CONFIG_PATHVLIST_EXTRA_ALLOC >=4 (guard values without bound checking)
#endif

static inline int clpPointInPolygon( const clpPoint64 *pt, int64_t *vertexlist, size_t vertexcount )
{
#if CLP_CONFIG_PATHVLIST_EXTRA_ALLOC >= 4
  return pwPointInContour64_Padding( vertexlist, vertexcount, pt->x, pt->y );
#else
  return pwPointInContour64( vertexlist, vertexcount, pt->x, pt->y );
#endif
}

static inline int clpPath1InsidePath2( const clpOutRec *or1, const clpOutRec *or2 )
{
  int result;
  clpOutPt *op;
  result = PW_RESULT_INSIDE;
  op = or1->pts;
  do
  {
    result = clpPointInPolygon( &op->pt, or2->path.vertexlist, or2->path.vertexcount );

#if 0
printf( "  Path1InsidePath2: %d\n", (int)result );
#endif

    if( result != PW_RESULT_ON )
      break;
    op = op->next;
  } while( op != or1->pts );

#if 0
printf( "Path1InsidePath2: %d\n", (int)result );
#endif

  if( result == PW_RESULT_ON )
    return ( fabs( clpGetOutPtArea( op, 0 ) ) < fabs( clpGetOutPtArea( or2->pts, 0 ) ) );
  else
    return ( result == PW_RESULT_INSIDE );
}


static inline void clpSetShapeTreeHoleFlags( clpTree *shapetree, int parentholeflag )
{
  int holeflag;
  size_t childindex;

  holeflag = 0;
  shapetree->holeflag = 0;
  if( shapetree->parent )
  {
    holeflag = parentholeflag ^ 1;
    shapetree->holeflag = holeflag ^ 1;
  }
  for( childindex = 0 ; childindex < shapetree->childcount ; childindex++ )
    clpSetShapeTreeHoleFlags( shapetree->childlist[ childindex ], holeflag );
  return;
}


////


static inline int clpBboxIsEmpty( int64_t *bbox )
{
  return ( bbox[0] >= bbox[1] );
}

static inline void clpGetContourBbox( int64_t *bbox, int64_t *vertexlist, size_t vertexcount )
{
  size_t vertexindex;
  int64_t vx, vy, bbminx, bbmaxx, bbminy, bbmaxy;
  bbminx = INT64_MAX;
  bbmaxx = INT64_MIN;
  bbminy = INT64_MAX;
  bbmaxy = INT64_MIN;
  /* TODO: Optimize! _mm_cmpgt_epi64() and _mm_blendv_ps() ~ no _mm_max_epi64() without AVX512 */
  for( vertexindex = 0 ; vertexindex < vertexcount ; vertexindex++ )
  {
    vx = vertexlist[ ( vertexindex * 2 ) + 0 ];
    vy = vertexlist[ ( vertexindex * 2 ) + 1 ];
    if( vx < bbminx )
      bbminx = vx;
    if( vx > bbmaxx )
      bbmaxx = vx;
    if( vy < bbminy )
      bbminy = vy;
    if( vy > bbmaxy )
      bbmaxy = vy;
  }
  bbox[0] = bbminx;
  bbox[1] = bbmaxx;
  bbox[2] = bbminy;
  bbox[3] = bbmaxy;
  return;
}

static inline int clpTestBboxInclusion( int64_t *parentbbox, int64_t *childbbox )
{
  /* TODO: Optimize! Individual branches are unpredictable, solve with SIMD to have just _one_ predictable branch */
  if( parentbbox[0] > childbbox[0] )
    return 0;
  if( parentbbox[1] < childbbox[1] )
    return 0;
  if( parentbbox[2] > childbbox[2] )
    return 0;
  if( parentbbox[3] < childbbox[3] )
    return 0;
  return 1;
}


static int clpDeepCheckOwner( clpClipper *clp, clpOutRec *outrec, clpOutRec *owner )
{
  size_t splitindex, splitcount;
  int insideownerbboxflag;
  clpOutRec *split;

  if( clpBboxIsEmpty( owner->bbox ) )
    clpGetContourBbox( owner->bbox, owner->path.vertexlist, owner->path.vertexcount );
  insideownerbboxflag = clpTestBboxInclusion( owner->bbox, outrec->bbox );

#if 0
printf( "DeepCheckOwner Bbox: %d %d %d %d\n", (int)owner->bbox[0], (int)owner->bbox[1], (int)owner->bbox[2], (int)owner->bbox[3] );
#endif

  /* While looking for the correct owner, check the owner's */
  /* splits **before** checking the owner itself because */
  /* splits can occur internally, and checking the owner */
  /* first would miss the inner split's true ownership */

  splitcount = clpOutRecVector_GetCount( &owner->splits );

  if( splitcount )
  { 
    for( splitindex = 0 ; splitindex < splitcount ; splitindex++ )
    {
      split = clpOutRecVector_Read( &owner->splits, splitindex );

#if 0
printf( "Test split %d\n", (int)split->idx );
#endif

      split = clpGetRealOutRec( split );
      if( !split || ( split->idx <= owner->idx ) || ( split == outrec ) )
        continue;
      if( clpOutRecVector_GetCount( &split->splits ) && clpDeepCheckOwner( clp, outrec, split ) )
        return 1;

      if( split->path.vertexcount == 0 )
        clpBuildPath( split->pts, clp->reversesolutionflag, 0, &split->path );
      if( clpBboxIsEmpty( split->bbox ) )
        clpGetContourBbox( split->bbox, split->path.vertexlist, split->path.vertexcount );

      if( clpTestBboxInclusion( split->bbox, outrec->bbox ) && clpPath1InsidePath2( outrec, split ) )
      {
        outrec->owner = split;
        return 1;
      }
    }
  }

#if 0
printf( "Foobar %d ~ Loop is_inside_owner_bounds %d\n", ( owner != outrec->owner ), (int)insideownerbboxflag );
#endif

  /* Only continue past here when not inside recursion */
  if( owner != outrec->owner )
    return 0;

  for( ; ; )
  {
    if( insideownerbboxflag && clpPath1InsidePath2( outrec, outrec->owner ) )
      return 1;
    /* Otherwise keep trying with owner's owner  */
    outrec->owner = outrec->owner->owner;
    if( !outrec->owner )
      return 1;
    insideownerbboxflag = clpTestBboxInclusion( (outrec->owner)->bbox, outrec->bbox );
  }
  return 1;
}


void clpBuildTree( clpClipper *clp, clpTree *shapetree, clpShape *shapeopen )
{
  size_t openindex, outrecindex, outreclistcount, childindex;
  clpOutRec *outrec, *tmp;
  size_t tmp_idx;
  clpTree *parentshapetree;
  clpTree *newshapetree;

  shapetree->path.vertexcount = 0;
  shapetree->path.vertexlist = 0;
  shapetree->parent = 0;
  shapetree->childcount = 0;
  shapetree->childlist = 0;
  shapetree->holeflag = 0;
  openindex = 0;
  shapeopen->pathcount = 0;
  shapeopen->pathlist = 0;
  outreclistcount = clpOutRecVector_GetCount( &clp->outreclist );
  if( clp->openpathsflag )
  {
    shapeopen->pathcount = 0;
    shapeopen->pathlist = (clpPath *)malloc( outreclistcount * sizeof(clpPath) );
  }

  for( outrecindex = 0 ; outrecindex < outreclistcount ; outrecindex++ )
  {
    outrec = clpOutRecVector_Read( &clp->outreclist, outrecindex );
    if( !outrec || !outrec->pts )
      continue;

    outrec->path.vertexcount = 0;
    outrec->path.vertexlist = 0;
    clpBboxSetEmpty( outrec->bbox );

    if( outrec->openflag )
    {
      if( clp->openpathsflag )
      {
        if( clpBuildPath( outrec->pts, clp->reversesolutionflag, 1, &shapeopen->pathlist[openindex] ) )
          openindex++;
      }
      continue;
    }

    if( !clpBuildPath( outrec->pts, clp->reversesolutionflag, 0, &outrec->path ) )
      continue;

    if( clpBboxIsEmpty( outrec->bbox ) )
      clpGetContourBbox( outrec->bbox, outrec->path.vertexlist, outrec->path.vertexcount );

#if 0
printf( "Bbox: %d %d %d %d\n", (int)outrec->bbox[0], (int)outrec->bbox[1], (int)outrec->bbox[2], (int)outrec->bbox[3] );
#endif

    outrec->owner = clpGetRealOutRec( outrec->owner );

#if 0
printf( "Outrec, owner %d, idx %d\n", (outrec->owner!=0), (int)outrec->idx );
#endif

    if( outrec->owner )
      clpDeepCheckOwner( clp, outrec, outrec->owner );

#if 0
printf( "Outrec, owner %d, idx %d\n", (outrec->owner!=0), (int)outrec->idx );
#endif

    /* Swap the order when a child preceeds its owner */
    /* (because owners must preceed children in polytrees) */
    if( outrec->owner && ( outrec->idx < (outrec->owner)->idx ) )
    {
      tmp = outrec->owner;
      clpOutRecVector_Write( &clp->outreclist, (outrec->owner)->idx, outrec );
      clpOutRecVector_Write( &clp->outreclist, outrec->idx, tmp );
      tmp_idx = outrec->idx;
      outrec->idx = tmp->idx;
      tmp->idx = tmp_idx;
      outrec = tmp;
      outrec->owner = clpGetRealOutRec( outrec->owner );
      /* Verify, is check required? */
      if( outrec->path.vertexcount == 0 )
        clpBuildPath( outrec->pts, clp->reversesolutionflag, 0, &outrec->path );
      if( clpBboxIsEmpty( outrec->bbox ) )
        clpGetContourBbox( outrec->bbox, outrec->path.vertexlist, outrec->path.vertexcount );
      if( outrec->owner )
        clpDeepCheckOwner( clp, outrec, outrec->owner );
    }

    if( outrec->owner && (outrec->owner)->shapetree )
      parentshapetree = (outrec->owner)->shapetree;
    else
      parentshapetree = shapetree;

#if 0
printf( "Outrec, owner %d, parent is root? %d\n", (outrec->owner!=0), ( parentshapetree == shapetree ) );
#endif

    childindex = parentshapetree->childcount++;
    newshapetree = (clpTree *)malloc( sizeof(clpTree) );
    parentshapetree->childlist = (clpTree **)realloc( parentshapetree->childlist, parentshapetree->childcount * sizeof(clpTree *) );
    parentshapetree->childlist[ childindex ] = newshapetree;
    newshapetree->path = outrec->path;
    newshapetree->parent = parentshapetree;
    newshapetree->childcount = 0;
    newshapetree->childlist = 0;
    newshapetree->holeflag = 0;

    outrec->shapetree = newshapetree;
  }

  clpSetShapeTreeHoleFlags( shapetree, 0 );

  return;
}


////


clpClipper *clpCreateClipper( int orientation_is_reversed )
{
  clpClipper *clp;
  clp = (clpClipper *)malloc( sizeof(clpClipper) );
  memset( (void *)clp, 0, sizeof(clpClipper) );

  clp->cliptype = CLP_CLIPTYPE_NONE;
  clp->fillrule = CLP_FILLRULE_EVENODD;
  clp->fillpos = CLP_FILLRULE_POSITIVE;
  clp->scanbeambottomy = 0;
  clp->openpathsflag = 0;
  clp->minimalistsorted = 0;
  clp->polytreeflag = 0;
  clp->succeessflag = 1;
  clp->activeedgelist = 0;
  clp->sortededgelist = 0;
  clp->horzjoinerlist = 0;
  clp->orientationreversedflag = orientation_is_reversed;
  if( clp->orientationreversedflag )
    clp->fillpos = CLP_FILLRULE_NEGATIVE;
  clp->preservecollinearflag = 1;
  clp->reversesolutionflag = 0;

  clpLocalMinimaVector_Init( &clp->minimalist, 0 );

  clpOutRecVector_Init( &clp->outreclist, 0 );

  clpIntersectNodeVector_Init( &clp->intersectnodelist, 0 );

  clpJoinerVector_Init( &clp->joinerlist, 0 );

  clpVertexList_Init( &clp->vertexlist, 0 );

  mmHeap64sbInit( &clp->scanlinelist, 0 );

#if CLP_CONFIG_ENABLE_SLAB_ALLOCATION
  mmSlabInit( &clp->slaboutpt, sizeof(clpOutPt), 4096, 0x10 );
  mmSlabInit( &clp->slablm, sizeof(clpLocalMinima), 1024, 0x10 );
  mmSlabInit( &clp->slaboutrec, sizeof(clpOutRec), 1024, 0x10 );
  mmSlabInit( &clp->slabjoiner, sizeof(clpJoiner), 512, 0x10 );
  mmSlabInit( &clp->slabactive, sizeof(clpActiveEdge), 4096, 0x10 );
#endif

  return clp;
}

void clpDestroyClipper( clpClipper *clp )
{
  clpClear( clp );
#if CLP_CONFIG_ENABLE_SLAB_ALLOCATION
  mmSlabFreeAll( &clp->slaboutpt );
  mmSlabFreeAll( &clp->slablm );
  mmSlabFreeAll( &clp->slaboutrec );
  mmSlabFreeAll( &clp->slabjoiner );
  mmSlabFreeAll( &clp->slabactive );
#endif
  free( clp );
  return;
}


void clpSetPreserveCollinear( clpClipper *clp, int preservecollinearflag )
{
  clp->preservecollinearflag = preservecollinearflag;
  return;
}

void clpSetReverseSolution( clpClipper *clp, int reversesolutionflag )
{
  clp->reversesolutionflag = reversesolutionflag;
  return;
}


////


void clpAddShape( clpClipper *clp, const clpShape *shape, uint8_t pathtype, uint8_t openflag )
{
  size_t pathindex, vertexindex, vertexcount;
  size_t total_vertex_count, validcount;
  clpVertex *vertices, *v;
  clpVertex *v0, *curv, *prevv;
  int upflag, upflag0;
  clpPath *path;
  clpPoint64 pt;

  if( openflag )
    clp->openpathsflag = 1;
  clp->minimalistsorted = 0;

  total_vertex_count = 0;
  for( pathindex = 0 ; pathindex < shape->pathcount ; pathindex++ )
    total_vertex_count += shape->pathlist[ pathindex ].vertexcount;
  if( total_vertex_count == 0 )
    return;

  vertices = (clpVertex *)malloc( total_vertex_count * sizeof(clpVertex) );
  clpVertexList_Push( &clp->vertexlist, vertices );

  v = vertices;
  vertexindex = 0;
  for( pathindex = 0 ; pathindex < shape->pathcount ; pathindex++ )
  {
    path = &shape->pathlist[ pathindex ];
    vertexcount = path->vertexcount;
    if( !vertexcount )
      continue;

    /* For each path create a circular double linked list of vertices */
    v0 = v;
    curv = v;
    prevv = 0;
    v->prev = 0;
    validcount = 0;
    for( vertexindex = 0 ; vertexindex < vertexcount ; vertexindex++ )
    {
      pt.x = path->vertexlist[ ( 2 * vertexindex ) + 0 ];
      pt.y = path->vertexlist[ ( 2 * vertexindex ) + 1 ];
      if( prevv )
      {
        if( clpPoint64Equal( &prevv->pt, &pt ) )
          continue;
        prevv->next = curv;
      }
      curv->prev = prevv;
      curv->pt = pt;
      curv->flags = 0;
      prevv = curv++;
      validcount++;
    }

    if( !prevv || !prevv->prev )
      continue;
    if( !openflag && clpPoint64Equal( &prevv->pt, &v0->pt ) )
      prevv = prevv->prev;
    prevv->next = v0;
    v0->prev = prevv;
    v = curv; /* Get ready for next path */

    if( ( validcount < 2 ) || ( ( validcount == 2 ) && !openflag ) )
      continue;

    /* Find and assign local minima */
    if( openflag )
    {
      curv = v0->next;
      while( ( curv != v0 ) && ( curv->pt.y == v0->pt.y ) )
        curv = curv->next;
      upflag = ( curv->pt.y <= v0->pt.y );
      if( upflag )
      {
        v0->flags = CLP_VERTEX_FLAG_OPENSTART;
        clpAddLocMin( clp, v0, pathtype, 1 );
      }
      else
        v0->flags = CLP_VERTEX_FLAG_OPENSTART | CLP_VERTEX_FLAG_LOCALMAX;
    }
    else
    {
      prevv = v0->prev;
      while( ( prevv != v0 ) && ( prevv->pt.y == v0->pt.y ) )
        prevv = prevv->prev;
      if( prevv == v0 )
        continue; /* Only open paths can be completely flat */
      upflag = ( prevv->pt.y > v0->pt.y );
    }

    upflag0 = upflag;
    prevv = v0;
    curv = v0->next;
    while( curv != v0 )
    {
      if( ( curv->pt.y > prevv->pt.y ) && upflag )
      {
        prevv->flags |= CLP_VERTEX_FLAG_LOCALMAX;
        upflag = 0;
      }
      else if( ( curv->pt.y < prevv->pt.y ) && !upflag )
      {
        upflag = 1;
        clpAddLocMin( clp, prevv, pathtype, openflag );
      }
      prevv = curv;
      curv = curv->next;
    }

    if( openflag )
    {
      prevv->flags |= CLP_VERTEX_FLAG_OPENEND;
      if( upflag )
        prevv->flags |= CLP_VERTEX_FLAG_LOCALMAX;
      else
        clpAddLocMin( clp, prevv, pathtype, openflag );
    }
    else if( upflag != upflag0 )
    {
      if( upflag0 )
        clpAddLocMin( clp, prevv, pathtype, 0 );
      else
        prevv->flags |= CLP_VERTEX_FLAG_LOCALMAX;
    }
  }

  return;
}

void clpAddPath( clpClipper *clp, int64_t *vertexlist, size_t vertexcount, uint8_t pathtype, uint8_t openflag )
{
  clpShape shape;
  clpPath path;
  shape.pathcount = 1;
  shape.pathlist = &path;
  path.vertexcount = vertexcount;
  path.vertexlist = vertexlist;
  clpAddShape( clp, &shape, pathtype, openflag );
  return;
}


////


int clpClipperExecute( clpClipper *clp, int cliptype, int fillrule, clpShape *shape )
{
  memset( (void *)shape, 0, sizeof(clpShape) );
  if( clpExecuteInternal( clp, cliptype, fillrule, 0 ) )
    clpBuildShape( clp, shape, 0 );
  clpCleanUp( clp );
  return clp->succeessflag;
}


int clpClipperExecute_RetOpen( clpClipper *clp, int cliptype, int fillrule, clpShape *shape, clpShape *shapeopen )
{
  memset( (void *)shape, 0, sizeof(clpShape) );
  memset( (void *)shapeopen, 0, sizeof(clpShape) );
  if( clpExecuteInternal( clp, cliptype, fillrule, 0 ) )
    clpBuildShape( clp, shape, shapeopen );
  clpCleanUp( clp );
  return clp->succeessflag;
}


int clpClipperExecute_RetTree( clpClipper *clp, int cliptype, int fillrule, clpTree *shapetree, clpShape *shapeopen )
{
  memset( (void *)shapetree, 0, sizeof(clpTree) );
  memset( (void *)shapeopen, 0, sizeof(clpShape) );
  if( clpExecuteInternal( clp, cliptype, fillrule, 1 ) )
    clpBuildTree( clp, shapetree, shapeopen );
  clpCleanUp( clp );
  return clp->succeessflag;
}


////


void clpFreeShape( clpShape *shape )
{
  size_t pathindex, pathcount;
  clpPath *path;

  pathcount = shape->pathcount;
  for( pathindex = 0 ; pathindex < pathcount ; pathindex++ )
  {
    path = &shape->pathlist[ pathindex ];
    free( path->vertexlist );
  }
  free( shape->pathlist );
  memset( shape, 0, sizeof(clpShape) );
  return;
}


static void clpFreeTreeChild( clpTree *tree )
{
  size_t childindex, childcount;
  childcount = tree->childcount;
  for( childindex = 0 ; childindex < childcount ; childindex++ )
    clpFreeTreeChild( tree->childlist[ childindex ] );
  free( tree->childlist );
  free( tree->path.vertexlist );
  free( tree );
  return;
}

void clpFreeTree( clpTree *tree )
{
  size_t childindex, childcount;
  childcount = tree->childcount;
  for( childindex = 0 ; childindex < childcount ; childindex++ )
    clpFreeTreeChild( tree->childlist[ childindex ] );
  free( tree->childlist );
  free( tree->path.vertexlist );
  memset( tree, 0, sizeof(clpTree) );
  return;
}


////


void clpConvDoubleToInt64( double *dstvl64, int64_t *srcvld, size_t vertexcount, double convfactor )
{
  size_t vertexindex;
  double convfactorinv;
  convfactorinv = 1.0;
  if( convfactor )
    convfactorinv = 1.0 / convfactor;
  for( vertexindex = 0 ; vertexindex < vertexcount ; vertexindex++ )
    dstvl64[ vertexindex ] = (int64_t)round( srcvld[ vertexindex ] * convfactorinv );
  return;
}

void clpConvInt64ToDouble( double *dstvld, int64_t *srcvl64, size_t vertexcount, double convfactor )
{
  size_t vertexindex;
  for( vertexindex = 0 ; vertexindex < vertexcount ; vertexindex++ )
    dstvld[ vertexindex ] = convfactor * (double)srcvl64[ vertexindex ];
  return;
}


////


int clpTestPointWithinShape( clpShape *shape, int64_t pointx, int64_t pointy )
{
  int parity;
  size_t pathindex;
  clpPath *path;
  parity = 0;
  for( pathindex = 0 ; pathindex < shape->pathcount ; pathindex++ )
  {
    path = &shape->pathlist[ pathindex ];
    if( pwPointInContour64( path->vertexlist, path->vertexcount, pointx, pointy ) != PW_RESULT_OUTSIDE )
      parity ^= 1;
  }
  return parity;
}