/*******************************************************************************
*                                                                              *
* 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 "polyclip.h"

#include "polyarea.h"


////


#ifndef M_PI
 #define M_PI (3.141592653589793238)
#endif

#define CLP_ARC_TOLERANCE (0.25)
#define CLP_EPSILON (1e-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;
}


#define MM_VECTOR_ITEMTYPE clpPoint64
#define MM_VECTOR_TYPENAME clpPointVector
#define MM_VECTOR_PREFIX clpPointVector_
#define MM_VECTOR_OPTION_ORDERED (1)
#define MM_VECTOR_OPTION_DEBUG_CALLS (0)
#include "mmvector.h"

#define MM_VECTOR_ITEMTYPE clpPath
#define MM_VECTOR_TYPENAME clpPathVector
#define MM_VECTOR_PREFIX clpPathVector_
#define MM_VECTOR_OPTION_ORDERED (1)
#define MM_VECTOR_OPTION_DEBUG_CALLS (0)
#include "mmvector.h"

#define MM_VECTOR_ITEMTYPE clpPointD
#define MM_VECTOR_TYPENAME clpPathD
#define MM_VECTOR_PREFIX clpPathD_
#define MM_VECTOR_OPTION_ORDERED (1)
#define MM_VECTOR_OPTION_DEBUG_CALLS (0)
#include "mmvector.h"


typedef struct
{
  /* Direct input shape from user */
  clpShape inputshape;
  /* Where we write our offseted vertices, before creating the final shape */
  clpPointVector midpath;
  /* Final shape for group */
  clpPathVector pathlist;
  /* Same pathlist as above, packed into shape struct */
  clpShape shape;

  uint8_t reversedflag;
  uint8_t jointype;
  uint8_t endtype;
  /* If non-zero, inputshape pathlist was allocated internally ~ it will be freed at the end */
  uint8_t inputpathlistallocflag;
} clpPathGroup;


#define MM_VECTOR_ITEMTYPE clpPathGroup
#define MM_VECTOR_TYPENAME clpPathGroupVector
#define MM_VECTOR_PREFIX clpPathGroupVector_
#define MM_VECTOR_OPTION_ORDERED (1)
#define MM_VECTOR_OPTION_DEBUG_CALLS (0)
#include "mmvector.h"


struct clpOffset
{
  double offsetdelta;
  double miterthreshold;
  double miterlimit;
  double roundradiansteps;
  clpPathD norms;

  clpPathVector pathlist;
  clpShape finalshape;

  clpPathGroupVector grouplist;

  // MergeGroups: A path group is one or more paths added via the AddPath or
  // AddPaths methods. By default these path groups will be offset
  // independently of other groups and this may cause overlaps (intersections).
  // However, when MergeGroups is enabled, any overlapping offsets will be
  // merged (via a clipping union operation) to remove overlaps.
  uint8_t jointype;
  uint8_t mergegroupsflag;
  uint8_t preservecollinearflag;
  uint8_t reversesolutionflag;
};


////


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


static inline size_t clpGetLowestPolygonIdx( clpShape *shape )
{
  size_t pathindex, vertexindex, vertexcount, bottompath;
  int64_t *vertexlist;
  clpPoint64 lp;

  bottompath = 0;
  lp = clpBuildPoint64( 0, INT64_MIN );
  for( pathindex = 0 ; pathindex < shape->pathcount ; pathindex++ )
  {
    vertexcount = shape->pathlist[pathindex].vertexcount;
    vertexlist = shape->pathlist[pathindex].vertexlist;
    for( vertexindex = 0 ; vertexindex < vertexcount ; vertexindex++ )
    {
      if( ( vertexlist[1] > lp.y ) || ( ( vertexlist[1] == lp.y ) && ( vertexlist[0] < lp.x ) ) )
      {
        bottompath = pathindex;
        lp = clpBuildPoint64( vertexlist[0], vertexlist[1] );
      }
      vertexlist += 2;
    }
  }
  return bottompath;
}



static inline clpPointD clpGetNormal( int64_t pt0x, int64_t pt0y, int64_t pt1x, int64_t pt1y )
{
  double dx, dy, inverse_hypot;
  dx = (double)( pt1x - pt0x );
  dy = (double)( pt1y - pt0y );
  inverse_hypot = 1.0 / hypot( dx, dy );
  dx *= inverse_hypot;
  dy *= inverse_hypot;
  return clpBuildPointD( dy, -dx );
}

static inline int clpIsClosedPath( int et )
{
  return ( ( et == CLP_ENDTYPE_POLYGON ) || ( et == CLP_ENDTYPE_JOINED ) );
}


////


static void clpBuildNormals( clpOffset *clpoffset, int64_t *vertexlist, size_t vertexcount )
{
  size_t vertexindex;
  int64_t ptx, pty, lastx, lasty;
  clpPointD pd;
  clpPathD_Resize( &clpoffset->norms, vertexcount, clpPathD_GetAlloc( &clpoffset->norms ) );
  clpPathD_Empty( &clpoffset->norms );
  if( vertexcount == 0 )
    return;
  lastx = vertexlist[ 0 ];
  lasty = vertexlist[ 1 ];
  for( vertexindex = 1 ; vertexindex < vertexcount ; vertexindex++ )
  {
    ptx = vertexlist[ ( vertexindex * 2 ) + 0 ];
    pty = vertexlist[ ( vertexindex * 2 ) + 1 ];
    pd = clpGetNormal( lastx, lasty, ptx, pty );
    clpPathD_Push( &clpoffset->norms, pd );
    lastx = ptx;
    lasty = pty;
  }
  ptx = vertexlist[ 0 ];
  pty = vertexlist[ 1 ];
  pd = clpGetNormal( lastx, lasty, ptx, pty );
  clpPathD_Push( &clpoffset->norms, pd );
  return;
}

static void clpDoSquare( clpOffset *clpoffset, clpPathGroup *group, int64_t *vertexlist, size_t j, size_t k )
{
  int64_t x, y;
  int64_t pathx, pathy;
  clpPointD pt0, pt1;

  pt0 = clpPathD_Read( &clpoffset->norms, j );
  pt1 = clpPathD_Read( &clpoffset->norms, k );

  pathx = vertexlist[ ( j * 2 ) + 0 ];
  pathy = vertexlist[ ( j * 2 ) + 1 ];

  if( clpoffset->offsetdelta > 0 )
  {
#if CLP_DEBUG_VERBOSE
    printf( "&& %lld + %.16f * ( %.16f - %.16f )\n", (long long)pathx, clpoffset->offsetdelta, pt1.x, pt1.y );
    printf( "&& %lld + %.16f * ( %.16f + %.16f )\n", (long long)pathy, clpoffset->offsetdelta, pt1.y, pt1.x );
#endif
    x = pathx + round( clpoffset->offsetdelta * ( pt1.x - pt1.y ) );
    y = pathy + round( clpoffset->offsetdelta * ( pt1.y + pt1.x ) );
    clpPointVector_Push( &group->midpath, clpBuildPoint64( x, y ) );
#if CLP_DEBUG_VERBOSE
    printf( "DoSquare ~ Push %lld %lld\n", (long long)x, (long long)y );
#endif
    x = pathx + round( clpoffset->offsetdelta * ( pt0.x + pt0.y ) );
    y = pathy + round( clpoffset->offsetdelta * ( pt0.y - pt0.x ) );
    clpPointVector_Push( &group->midpath, clpBuildPoint64( x, y ) );
#if CLP_DEBUG_VERBOSE
    printf( "DoSquare ~ Push %lld %lld\n", (long long)x, (long long)y );
#endif
  }
  else
  {
#if CLP_DEBUG_VERBOSE
    printf( "&& %lld + %.16f * ( %.16f + %.16f )\n", (long long)pathx, clpoffset->offsetdelta, pt1.x, pt1.y );
    printf( "&& %lld + %.16f * ( %.16f - %.16f )\n", (long long)pathy, clpoffset->offsetdelta, pt1.y, pt1.x );
#endif
    x = pathx + round( clpoffset->offsetdelta * ( pt1.x + pt1.y ) );
    y = pathy + round( clpoffset->offsetdelta * ( pt1.y - pt1.x ) );
    clpPointVector_Push( &group->midpath, clpBuildPoint64( x, y ) );
#if CLP_DEBUG_VERBOSE
    printf( "DoSquare ~ Push %lld %lld\n", (long long)x, (long long)y );
#endif
    x = pathx + round( clpoffset->offsetdelta * ( pt0.x - pt0.y ) );
    y = pathy + round( clpoffset->offsetdelta * ( pt0.y + pt0.x ) );
    clpPointVector_Push( &group->midpath, clpBuildPoint64( x, y ) );
#if CLP_DEBUG_VERBOSE
    printf( "DoSquare ~ Push %lld %lld\n", (long long)x, (long long)y );
#endif
  }

  return;
}

static void clpDoMiter( clpOffset *clpoffset, clpPathGroup *group, int64_t *vertexlist, size_t j, size_t k, double cos_a )
{
  int64_t x, y;
  int64_t pathx, pathy;
  double q;
  clpPointD pt0, pt1;

  pt0 = clpPathD_Read( &clpoffset->norms, j );
  pt1 = clpPathD_Read( &clpoffset->norms, k );

  pathx = vertexlist[ ( j * 2 ) + 0 ];
  pathy = vertexlist[ ( j * 2 ) + 1 ];

  q = clpoffset->offsetdelta / ( cos_a + 1.0 );
  x = pathx + round( ( pt1.x + pt0.x ) * q );
  y = pathy + round( ( pt1.y + pt0.y ) * q );
  clpPointVector_Push( &group->midpath, clpBuildPoint64( x, y ) );
#if CLP_DEBUG_VERBOSE
  printf( "DoMiter ~ Push %lld %lld\n", (long long)x, (long long)y );
#endif
  return;
}

static void clpDoRound( clpOffset *clpoffset, clpPathGroup *group, const clpPoint64 *pt, const clpPointD *norm1, const clpPointD *norm2, double angle, int minstepcount )
{
  int i;
  int64_t x, y;
  clpPointD pt2;
  int steps;
  double step_sin, step_cos;

  /* Even though angle may be negative this is a convex join */
  pt2 = clpBuildPointD( norm2->x * clpoffset->offsetdelta, norm2->y * clpoffset->offsetdelta );
  steps = (int)( round( clpoffset->roundradiansteps * fabs( angle ) + 0.501 ) );
  if( steps < minstepcount )
    steps = minstepcount;
  x = pt->x + round( pt2.x );
  y = pt->y + round( pt2.y );
  clpPointVector_Push( &group->midpath, clpBuildPoint64( x, y ) );
#if CLP_DEBUG_VERBOSE
  printf( "DoRound ~ Push %lld %lld\n", (long long)x, (long long)y );
#endif

  step_sin = sin( angle / steps );
  step_cos = cos( angle / steps );
  for( i = 1 ; i < steps ; i++ )
  {
    pt2 = clpBuildPointD( ( pt2.x * step_cos ) - ( step_sin * pt2.y ), ( pt2.x * step_sin ) + ( pt2.y * step_cos ) );
    x = pt->x + round( pt2.x );
    y = pt->y + round( pt2.y );
    clpPointVector_Push( &group->midpath, clpBuildPoint64( x, y ) );
#if CLP_DEBUG_VERBOSE
    printf( "DoRound ~ Push %lld %lld\n", (long long)x, (long long)y );
#endif
  }
  pt2 = clpBuildPointD( norm1->x * clpoffset->offsetdelta, norm1->y * clpoffset->offsetdelta );
  x = pt->x + round( pt2.x );
  y = pt->y + round( pt2.y );
  clpPointVector_Push( &group->midpath, clpBuildPoint64( x, y ) );
#if CLP_DEBUG_VERBOSE
  printf( "DoRound ~ Push %lld %lld\n", (long long)x, (long long)y );
#endif

  return;
}

static void clpOffsetPoint( clpOffset *clpoffset, clpPathGroup *group, int64_t *vertexlist, size_t j, size_t k )
{
  // A: angle between adjoining edges (on left side WRT winding direction).
  // A == 0 deg (or A == 360 deg): collinear edges heading in same direction
  // A == 180 deg: collinear edges heading in opposite directions (i.e. a 'spike')
  // sin(A) < 0: convex on left.
  // cos(A) > 0: angles on both left and right sides > 90 degrees
  int64_t x, y;
  int64_t pathx, pathy;
  int tinyangleflag;
  double sin_a, cos_a;
  clpPointD pt0, pt1;
  clpPoint64 p1, p2;

  pt0 = clpPathD_Read( &clpoffset->norms, j );
  pt1 = clpPathD_Read( &clpoffset->norms, k );

  sin_a = ( pt1.x * pt0.y ) - ( pt0.x * pt1.y );
  sin_a = fmax( fmin( sin_a, 1.0 ), -1.0 );
  cos_a = ( pt0.x * pt1.x ) + ( pt0.y * pt1.y );
  tinyangleflag = ( fabs( sin_a ) < 0.001 ) && ( cos_a > 0.0 );

  pathx = vertexlist[ ( j * 2 ) + 0 ];
  pathy = vertexlist[ ( j * 2 ) + 1 ];

  /* Almost no angle of deviation or it's concave */
  if( ( tinyangleflag ) || ( ( sin_a * clpoffset->offsetdelta ) < 0.0 ) )
  {
    x = pathx + round( pt1.x * clpoffset->offsetdelta );
    y = pathy + round( pt1.y * clpoffset->offsetdelta );
    p1 = clpBuildPoint64( x, y );
    x = pathx + round( pt0.x * clpoffset->offsetdelta );
    y = pathy + round( pt0.y * clpoffset->offsetdelta );
    p2 = clpBuildPoint64( x, y );
    clpPointVector_Push( &group->midpath, p1 );
#if CLP_DEBUG_VERBOSE
    printf( "OffsetPoint ~ Push %lld %lld\n", (long long)p1.x, (long long)p1.y );
#endif
    if( !clpPoint64Equal( &p1, &p2 ) )
    {
      clpPointVector_Push( &group->midpath, clpBuildPoint64( pathx, pathy ) );
#if CLP_DEBUG_VERBOSE
      printf( "OffsetPoint ~ Push %lld %lld\n", (long long)pathx, (long long)pathy );
#endif
      clpPointVector_Push( &group->midpath, p2 );
#if CLP_DEBUG_VERBOSE
      printf( "OffsetPoint ~ Push %lld %lld\n", (long long)p2.x, (long long)p2.y );
#endif
    }
  }
  else
  {
    switch( clpoffset->jointype )
    {
      case CLP_JOINTYPE_MITER:
        if( ( 1.0 + cos_a ) < clpoffset->miterthreshold )
          clpDoSquare( clpoffset, group, vertexlist, j, k );
        else
          clpDoMiter( clpoffset, group, vertexlist, j, k, cos_a );
        break;
      case CLP_JOINTYPE_SQUARE:
        if( cos_a >= 0.0 )
          clpDoMiter( clpoffset, group, vertexlist, j, k, cos_a );
        else
          clpDoSquare( clpoffset, group, vertexlist, j, k );
        break;
      default:
        p1 = clpBuildPoint64( pathx, pathy );
        clpDoRound( clpoffset, group, &p1, &pt0, &pt1, atan2( sin_a, cos_a ), 0 );
        break;
    }
  }
  return;
}

static void clpOffsetPolygon( clpOffset *clpoffset, clpPathGroup *group, int64_t *vertexlist, size_t vertexcount )
{
  size_t i, j;
  clpPath newpath;

  clpPointVector_Clear( &group->midpath );
  j = vertexcount - 1;
  for( i = 0 ; i < vertexcount ; i++ )
  {
    clpOffsetPoint( clpoffset, group, vertexlist, i, j );
    j = i;
  }
  newpath.vertexcount = clpPointVector_GetCount( &group->midpath );
  newpath.vertexlist = (int64_t *)clpPointVector_Steal( &group->midpath );
  clpPathVector_Push( &group->pathlist, newpath );

  return;
}

static inline void clpOffsetOpenJoined( clpOffset *clpoffset, clpPathGroup *group, int64_t *vertexlist, size_t vertexcount )
{
  size_t vertexindex, vertexswapstop, swapindex, vertexlast;
  int64_t swapx, swapy;

  clpOffsetPolygon( clpoffset, group, vertexlist, vertexcount );

  /* Revert the order of vertices, in a clunky way */
  vertexswapstop = vertexcount >> 1;
  vertexlast = vertexcount - 1;
  for( vertexindex = 0 ; vertexindex < vertexswapstop ; vertexindex++ )
  {
    swapx = vertexlist[ ( vertexindex * 2 ) + 0 ];
    swapy = vertexlist[ ( vertexindex * 2 ) + 1 ];
    swapindex = vertexlast - vertexindex;
    vertexlist[ ( vertexindex * 2 ) + 0 ] = vertexlist[ ( swapindex * 2 ) + 0 ];
    vertexlist[ ( vertexindex * 2 ) + 1 ] = vertexlist[ ( swapindex * 2 ) + 1 ];
    vertexlist[ ( swapindex * 2 ) + 0 ] = swapx;
    vertexlist[ ( swapindex * 2 ) + 1 ] = swapy;
  }

  clpBuildNormals( clpoffset, vertexlist, vertexcount );
  clpOffsetPolygon( clpoffset, group, vertexlist, vertexcount );
  return;
}

static inline void clpOffsetOpenPath( clpOffset *clpoffset, clpPathGroup *group, int64_t *vertexlist, size_t vertexcount, int end_type )
{
  size_t i, j, k;
  int64_t x, y;
  int64_t pathx, pathy;
  clpPointD normpt0, normpt1;
  clpPath newpath;
  clpPoint64 endpt;

  clpPointVector_Clear( &group->midpath );

  j = 0;
  for( i = 1 ; i < vertexcount - 1 ; i++ )
  {
    clpOffsetPoint( clpoffset, group, vertexlist, i, j );
    j = i;
  }

  j = clpPathD_GetCount( &clpoffset->norms ) - 1;
  k = j - 1;
  normpt0 = clpPathD_Read( &clpoffset->norms, k );
  normpt1 = clpBuildPointD( -normpt0.x, -normpt0.y );
  clpPathD_Write( &clpoffset->norms, j, normpt1 );

  pathx = vertexlist[ ( j * 2 ) + 0 ];
  pathy = vertexlist[ ( j * 2 ) + 1 ];

  switch( end_type )
  {
    case CLP_ENDTYPE_BUTT:
      x = pathx + round( normpt0.x * clpoffset->offsetdelta );
      y = pathy + round( normpt0.y * clpoffset->offsetdelta );
      clpPointVector_Push( &group->midpath, clpBuildPoint64( x, y ) );
#if CLP_DEBUG_VERBOSE
      printf( "OffsetOpenPath ~ Push %lld %lld\n", (long long)x, (long long)y );
#endif
      x = pathx - round( normpt0.x * clpoffset->offsetdelta );
      y = pathy - round( normpt0.y * clpoffset->offsetdelta );
      clpPointVector_Push( &group->midpath, clpBuildPoint64( x, y ) );
#if CLP_DEBUG_VERBOSE
      printf( "OffsetOpenPath ~ Push %lld %lld\n", (long long)x, (long long)y );
#endif
      break;
    case CLP_ENDTYPE_ROUND:
      endpt = clpBuildPoint64( pathx, pathy );
      clpDoRound( clpoffset, group, &endpt, &normpt1, &normpt0, M_PI, 0 );
      break;
    default:
      clpDoSquare( clpoffset, group, vertexlist, j, k );
      break;
  }

  /* Reverse normals, an open path has to be expanded on both sides */
  for( i = k ; i > 0 ; i-- )
  {
    normpt0 = clpPathD_Read( &clpoffset->norms, i-1 );
    clpPathD_Write( &clpoffset->norms, i, clpBuildPointD( -normpt0.x, -normpt0.y ) );
  }
  normpt0 = clpPathD_Read( &clpoffset->norms, 1 );
  clpPathD_Write( &clpoffset->norms, 0, clpBuildPointD( -normpt0.x, -normpt0.y ) );

  for( i = k ; i > 0 ; i-- )
  {
    clpOffsetPoint( clpoffset, group, vertexlist, i, j );
    j = i;
  }

  pathx = vertexlist[ 0 ];
  pathy = vertexlist[ 1 ];

  /* Cap the end point */
  switch( end_type )
  {
    case CLP_ENDTYPE_BUTT:
      normpt0 = clpPathD_Read( &clpoffset->norms, 1 );
      x = pathx + round( normpt0.x * clpoffset->offsetdelta );
      y = pathy + round( normpt0.y * clpoffset->offsetdelta );
      clpPointVector_Push( &group->midpath, clpBuildPoint64( x, y ) );
#if CLP_DEBUG_VERBOSE
      printf( "OffsetOpenPath ~ Push %lld %lld\n", (long long)x, (long long)y );
#endif
      x = pathx - round( normpt0.x * clpoffset->offsetdelta );
      y = pathy - round( normpt0.y * clpoffset->offsetdelta );
      clpPointVector_Push( &group->midpath, clpBuildPoint64( x, y ) );
#if CLP_DEBUG_VERBOSE
      printf( "OffsetOpenPath ~ Push %lld %lld\n", (long long)x, (long long)y );
#endif
      break;
    case CLP_ENDTYPE_ROUND:
      endpt = clpBuildPoint64( pathx, pathy );
      normpt0 = clpPathD_Read( &clpoffset->norms, 0 );
      normpt1 = clpPathD_Read( &clpoffset->norms, 1 );
      clpDoRound( clpoffset, group, &endpt, &normpt0, &normpt1, M_PI, 0 );
      break;
    default:
      clpDoSquare( clpoffset, group, vertexlist, 0, 1 );
      break;
  }

  newpath.vertexcount = clpPointVector_GetCount( &group->midpath );
  newpath.vertexlist = (int64_t *)clpPointVector_Steal( &group->midpath );
  clpPathVector_Push( &group->pathlist, newpath );
  return;
}

static inline double clpClosedPathArea( int64_t *vertexlist, size_t vertexcount, int orientation_is_reversed )
{
#if 1
  double area;
  area = paAccurateArea64_valid62( vertexlist, vertexcount );
#else
  size_t vertexindex, vertexlast, vertexstop;
  double area;
  int64_t v0x, v0y, v1x, v1y;

  if( vertexcount < 3 )
    return 0.0;
  area = 0.0;
  vertexlast = vertexcount - 1;
  v1x = vertexlist[ ( vertexlast * 2 ) + 0 ];
  v1y = vertexlist[ ( vertexlast * 2 ) + 1 ];
  vertexstop = vertexcount & ~1;
  for( vertexindex = 0 ; vertexindex < vertexstop ; vertexindex += 2 )
  {
    v0x = vertexlist[0];
    v0y = vertexlist[1];
    area += (double)( v1y + v0y ) * ( v1x - v0x );
    v1x = vertexlist[2];
    v1y = vertexlist[3];
    area += (double)( v0y + v1y ) * ( v0x - v1x );
    vertexlist += 4;
  }
  if( vertexcount & 1 )
  {
    v0x = vertexlist[0];
    v0y = vertexlist[1];
    area += (double)( v1y + v0y ) * ( v1x - v0x );
  }
  area *= 0.5;
#endif
  if( orientation_is_reversed )
    area = -area;
  return area;
}


static inline int64_t *clpCopyStripDuplicates( int64_t *vertexlist, size_t *inoutvertexcount, int closedpathflag )
{
  size_t vertexstart, vertexindex, vertexcount;
  int64_t lastx, lasty;
  int64_t *newvertexlist, *writevl;

  vertexcount = *inoutvertexcount;
  if( !vertexcount )
    return 0;
  newvertexlist = (int64_t *)malloc( vertexcount * 2 * sizeof(int64_t) );

  writevl = newvertexlist;
  if( closedpathflag )
  {
    lastx = vertexlist[ ( vertexcount * 2 ) - 2 ];
    lasty = vertexlist[ ( vertexcount * 2 ) - 1 ];
    vertexstart = 0;
  }
  else
  {
    writevl[0] = vertexlist[0];
    writevl[1] = vertexlist[1];
    lastx = vertexlist[0];
    lasty = vertexlist[1];
    writevl += 2;
    vertexlist += 2;
    vertexstart = 1;
  }

  for( vertexindex = vertexstart ; vertexindex < vertexcount ; vertexindex++, vertexlist += 2 )
  {
    if( ( vertexlist[0] != lastx ) || ( vertexlist[1] != lasty ) )
    {
      writevl[0] = vertexlist[0];
      writevl[1] = vertexlist[1];
      lastx = vertexlist[0];
      lasty = vertexlist[1];
      writevl += 2;
    }
  }

  vertexcount = ( writevl - newvertexlist ) >> 1;
  *inoutvertexcount = vertexcount;

  return newvertexlist;
}


static inline void clpDoGroupOffset( clpOffset *clpoffset, clpPathGroup *group, double delta )
{
  int closedpathflag;
  double area;
  size_t bottompath;
  clpPath *path;
  clpShape *inputshape;
  int64_t *vertexlist;
  size_t vertexcount;
  size_t pathindex, pathcount;
  clpShape *shape;
  clpClipper *clp;

  inputshape = &group->inputshape;

  if( group->endtype != CLP_ENDTYPE_POLYGON )
    delta = 0.5 * fabs( delta );
  closedpathflag = clpIsClosedPath( group->endtype );

  if( closedpathflag )
  {
    // the lowermost polygon must be an outer polygon. So we can use that as the
    // designated orientation for outer polygons (needed for tidy-up clipping)
    bottompath = clpGetLowestPolygonIdx( inputshape );
    // nb: don't use the default orientation here ...
    path = &inputshape->pathlist[ bottompath ];

    area = clpClosedPathArea( path->vertexlist, path->vertexcount, 0 );
#if CLP_DEBUG_VERBOSE
    printf( "%%%% AREA %f\n", area );
#endif
    if( area == 0.0 )
      return;
    group->reversedflag = ( area < 0.0 );
    if( group->reversedflag )
      delta = -delta;
  }
  else
    group->reversedflag = 0;

  clpoffset->offsetdelta = delta;
  clpoffset->jointype = group->jointype;

  closedpathflag = clpIsClosedPath( group->endtype );

  pathcount = inputshape->pathcount;

  for( pathindex = 0 ; pathindex < pathcount ; pathindex++ )
  {
    path = &inputshape->pathlist[ pathindex ];
    vertexcount = path->vertexcount;
    if( vertexcount == 0 )
      continue;
    vertexlist = clpCopyStripDuplicates( path->vertexlist, &vertexcount, closedpathflag );
    if( vertexcount == 1 ) /* single point - do we care about that case? */
    {
      clpPointD ptnorm;
      clpPoint64 ptvertex;
      clpPath newpath;
      clpPointVector_Clear( &group->midpath );
      ptnorm.x = 1.0;
      ptnorm.y = 0.0;
      ptvertex = clpBuildPoint64( vertexlist[0], vertexlist[1] );
      clpDoRound( clpoffset, group, &ptvertex, &ptnorm, &ptnorm, 2.0 * M_PI, 8 );
      newpath.vertexcount = clpPointVector_GetCount( &group->midpath ) - 1;
      newpath.vertexlist = (int64_t *)clpPointVector_Steal( &group->midpath );
      clpPathVector_Push( &group->pathlist, newpath );
    }
    else
    {
      clpBuildNormals( clpoffset, vertexlist, vertexcount );
      if( group->endtype == CLP_ENDTYPE_POLYGON )
        clpOffsetPolygon( clpoffset, group, vertexlist, vertexcount );
      else if( group->endtype == CLP_ENDTYPE_JOINED )
        clpOffsetOpenJoined( clpoffset, group, vertexlist, vertexcount );
      else
        clpOffsetOpenPath( clpoffset, group, vertexlist, vertexcount, group->endtype );
    }
    free( vertexlist );
  }

  group->shape.pathcount = clpPathVector_GetCount( &group->pathlist );
  group->shape.pathlist = clpPathVector_Steal( &group->pathlist );

  if( !clpoffset->mergegroupsflag )
  {
    /* Clean up self-intersections */
    clp = clpCreateClipper( 0 );

    clpSetPreserveCollinear( clp, 0 );
    /* The solution should retain the orientation of the input */
    clpSetReverseSolution( clp, ( clpoffset->reversesolutionflag != group->reversedflag ) );

    clpAddShape( clp, &group->shape, CLP_PATHTYPE_SUBJECT, 0 );

    if( group->reversedflag )
      clpClipperExecute( clp, CLP_CLIPTYPE_UNION, CLP_FILLRULE_NEGATIVE, &group->shape );
    else
      clpClipperExecute( clp, CLP_CLIPTYPE_UNION, CLP_FILLRULE_POSITIVE, &group->shape );

    clpDestroyClipper( clp );
  }

  shape = &group->shape;
  for( pathindex = 0 ; pathindex < shape->pathcount ; pathindex++ )
    clpPathVector_Push( &clpoffset->pathlist, shape->pathlist[ pathindex ] );

  free( group->shape.pathlist );
  group->shape.pathlist = 0;

  return;
}


////


void clpOffsetAddShape( clpOffset *clpoffset, clpShape *shape, int jointype, int endtype )
{
  clpPathGroup *pg;
  if( shape->pathcount == 0 )
    return;
  pg = clpPathGroupVector_Access( &clpoffset->grouplist, clpPathGroupVector_Add( &clpoffset->grouplist ) );
  memset( (void *)pg, 0, sizeof(clpPathGroup) );
  pg->reversedflag = 0;
  pg->jointype = jointype;
  pg->endtype = endtype;
  pg->inputshape = *shape;
  return;
}

void clpOffsetAddPath( clpOffset *clpoffset, int64_t *vertexlist, size_t vertexcount, int jointype, int endtype )
{
  clpPathGroup *pg;
  clpShape shape;
  clpPath *path;
  shape.pathcount = 1;
  shape.pathlist = (clpPath *)malloc( sizeof(clpPath) );
  path = shape.pathlist;
  path->vertexcount = vertexcount;
  path->vertexlist = vertexlist;
  clpOffsetAddShape( clpoffset, &shape, jointype, endtype );
  pg = clpPathGroupVector_Access( &clpoffset->grouplist, clpPathGroupVector_GetCount( &clpoffset->grouplist ) - 1 );
  pg->inputpathlistallocflag = 1;
  return;
}


////


#define CLP_OFFSET_DEBUG_STDOUT (0)
#define CLP_OFFSET_DEBUG_IMAGES (0)

#define CLP_OFFSET_DEBUG_IMAGE_SIZEX (4096)
#define CLP_OFFSET_DEBUG_IMAGE_SIZEY (4096)

#if CLP_OFFSET_DEBUG_IMAGES
 #include "img.h"
 #include "imgpng.h"
 #include "drawgraph.h"
#endif


#if CLP_OFFSET_DEBUG_STDOUT

static inline void clpDebugShape( clpShape *shape, char *shapename )
{
  size_t pathindex, vertexindex;
  double area;
  clpPath *clppath;
  printf( ">> Debug shape \"%s\", %d paths\n", shapename, (int)shape->pathcount );
  for( pathindex = 0 ; pathindex < shape->pathcount ; pathindex++ )
  {
    clppath = &shape->pathlist[ pathindex ];
    area = clpClosedPathArea( clppath->vertexlist, clppath->vertexcount, 0 );
    printf( "  Path %d ~ area %e\n", (int)pathindex, area );
    for( vertexindex = 0 ; vertexindex < clppath->vertexcount ; vertexindex++ )
      printf( "    Vertex %d : %lld %lld\n", (int)vertexindex, (long long)clppath->vertexlist[ 2*vertexindex + 0 ], (long long)clppath->vertexlist[ 2*vertexindex + 1 ] );
  }
  printf( "<< Debug shape \"%s\"\n", shapename );
  return;
}

#endif

#if CLP_OFFSET_DEBUG_IMAGES

typedef struct
{
  imgImage image;
  dgDmap *dmap;
  /* draw = ( input - base ) * scale */
  int64_t basex, basey;
  double scalex, scaley;
} clpOffsetDebugImage;

static inline void clpDebugShapeUpdateBbox( clpShape *shape, int64_t *bbox )
{
  size_t pathindex, vertexindex;
  int64_t px, py;
  clpPath *clppath;
  for( pathindex = 0 ; pathindex < shape->pathcount ; pathindex++ )
  {
    clppath = &shape->pathlist[ pathindex ];
    for( vertexindex = 0 ; vertexindex < clppath->vertexcount ; vertexindex++ )
    {
      px = clppath->vertexlist[ 2*vertexindex + 0 ];
      py = clppath->vertexlist[ 2*vertexindex + 1 ];
      if( px < bbox[0] )
        bbox[0] = px;
      if( px > bbox[1] )
        bbox[1] = px;
      if( py < bbox[2] )
        bbox[2] = py;
      if( py > bbox[3] )
        bbox[3] = py;
    }
  }
  return;
}

static inline void clpDebugInitTransform( clpOffsetDebugImage *debugimage, int64_t bbminx, int64_t bbmaxx, int64_t bbminy, int64_t bbmaxy, int drawsizex, int drawsizey, double extension )
{
  double scale;
  bbminx -= extension * (double)( bbmaxx - bbminx );
  bbmaxx += extension * (double)( bbmaxx - bbminx );
  bbminy -= extension * (double)( bbmaxy - bbminy );
  bbmaxy += extension * (double)( bbmaxy - bbminy );
  debugimage->basex = bbminx;
  debugimage->basey = bbminy;
  scale = fmin( (double)drawsizex / ( bbmaxx - bbminx ), (double)drawsizey / ( bbmaxy - bbminy ) );
  debugimage->scalex = scale;
  debugimage->scaley = scale;
  return;
}

static inline void clpDebugTransform2( float *out, clpOffsetDebugImage *debugimage, int64_t *v )
{
  out[0] = ( v[0] - debugimage->basex ) * debugimage->scalex;
  out[1] = ( v[1] - debugimage->basey ) * debugimage->scaley;
  out[2] = ( v[2] - debugimage->basex ) * debugimage->scalex;
  out[3] = ( v[3] - debugimage->basey ) * debugimage->scaley;
  return;
}

static void clpDebugDrawShape( clpOffsetDebugImage *debugimage, clpShape *shape, float *color )
{
  size_t pathindex, vertexindex;
  int64_t v[4];
  float line[4];
  clpPath *clppath;
  dgClearDmap( debugimage->dmap, 2.0f );
  for( pathindex = 0 ; pathindex < shape->pathcount ; pathindex++ )
  {
    clppath = &shape->pathlist[ pathindex ];
    if( !clppath->vertexcount )
      continue;
    for( vertexindex = 0 ; vertexindex < clppath->vertexcount ; vertexindex++ )
    {
      v[0] = clppath->vertexlist[ 2*vertexindex + 0 ];
      v[1] = clppath->vertexlist[ 2*vertexindex + 1 ];
      if( vertexindex < clppath->vertexcount-1 )
      {
        v[2] = clppath->vertexlist[ 2*vertexindex + 2 ];
        v[3] = clppath->vertexlist[ 2*vertexindex + 3 ];
      }
      else
      {
        v[2] = clppath->vertexlist[ 0 ];
        v[3] = clppath->vertexlist[ 1 ];
      }
      clpDebugTransform2( line, debugimage, v );
      dgAccumSegment( debugimage->dmap, line[0], line[1], line[2], line[3], 1.0f, 1.0f );
    }
  }
  dgRenderBlend( debugimage->image.data, debugimage->image.format.bytesperline, debugimage->dmap, color, 0.25f );
  return;
}


static inline void clpDebugInitImage( clpOffsetDebugImage *debugimage, int64_t *bbox, int drawsizex, int drawsizey, float extension )
{
  char *imagedata;
  size_t shapeindex, imagedataindex, imagedatasize;
  clpDebugInitTransform( debugimage, (double)bbox[0], (double)bbox[1], (double)bbox[2], (double)bbox[3], drawsizex, drawsizey, extension );
  imgNewImage( &debugimage->image, drawsizex, drawsizey, IMG_FORMAT_TYPE_RGBA8, 4 );
  imagedatasize = drawsizex * drawsizey * 4;
  imagedata = debugimage->image.data;
  for( imagedataindex = 0 ; imagedataindex < imagedatasize ; imagedataindex += 4 )
  {
    imagedata[imagedataindex+0] = 0;
    imagedata[imagedataindex+1] = 0;
    imagedata[imagedataindex+2] = 0;
    imagedata[imagedataindex+3] = 255;
  }
  debugimage->dmap = dgInitDmap( drawsizex, drawsizey );
  return;
}

static inline void clpDebugSaveImage( clpOffsetDebugImage *debugimage, char *path )
{
#if CLP_OFFSET_DEBUG_STDOUT
  printf( "!! Debug shape saved as \"%s\"\n", path );
#endif
  imgWritePngFile( path, &debugimage->image, 0.5 );
  return;
}

static inline void clpDebugFreeImage( clpOffsetDebugImage *debugimage )
{
  imgFree( &debugimage->image );
  dgFreeDmap( debugimage->dmap );
  return;
}

#endif


////



clpShape clpOffsetExecute( clpOffset *clpoffset, double delta )
{
  clpClipper *clp;
  int reversedflag;
  size_t groupindex, groupcount;
  clpShape workshape;
#if CLP_OFFSET_DEBUG_IMAGES
  char debugimagepath[256];
  static int debugimageindex = 0;
#endif

  clpPathVector_Clear( &clpoffset->pathlist );

  clpoffset->finalshape.pathcount = 0;
  clpoffset->finalshape.pathlist = 0;

#if CLP_OFFSET_DEBUG_STDOUT
  groupcount = clpPathGroupVector_GetCount( &clpoffset->grouplist );
  for( groupindex = 0 ; groupindex < groupcount ; groupindex++ )
    clpDebugShape( &(clpPathGroupVector_Access( &clpoffset->grouplist, groupindex )->inputshape), "Input" );
#endif

#if CLP_OFFSET_DEBUG_IMAGES
  int64_t bbox[4];
  clpOffsetDebugImage debugimage;
  bbox[0] = ((int64_t)1)<<62;
  bbox[1] = -((int64_t)1)<<62;
  bbox[2] = ((int64_t)1)<<62;
  bbox[3] = -((int64_t)1)<<62;
  for( groupindex = 0 ; groupindex < groupcount ; groupindex++ )
    clpDebugShapeUpdateBbox( &(clpPathGroupVector_Access( &clpoffset->grouplist, groupindex )->inputshape), bbox );
  bbox[0] -= delta;
  bbox[1] += delta;
  bbox[2] -= delta;
  bbox[3] += delta;
  clpDebugInitImage( &debugimage, bbox, CLP_OFFSET_DEBUG_IMAGE_SIZEX, CLP_OFFSET_DEBUG_IMAGE_SIZEY, 0.5 );
  for( groupindex = 0 ; groupindex < groupcount ; groupindex++ )
    clpDebugDrawShape( &debugimage, &(clpPathGroupVector_Access( &clpoffset->grouplist, groupindex )->inputshape), (float[]){ 128.0f, 0.0f, 0.0f, 255.0f } );
  snprintf( debugimagepath, 256, "debug-polyoffset/%05da.png", debugimageindex );
  clpDebugSaveImage( &debugimage, debugimagepath );
#endif

  clpoffset->miterthreshold = ( ( clpoffset->miterlimit <= 1 ) ? 2.0 : 2.0 / ( clpoffset->miterlimit * clpoffset->miterlimit ) );

  groupcount = clpPathGroupVector_GetCount( &clpoffset->grouplist );
  for( groupindex = 0 ; groupindex < groupcount ; groupindex++ )
    clpDoGroupOffset( clpoffset, clpPathGroupVector_Access( &clpoffset->grouplist, groupindex ), delta );

#if CLP_OFFSET_DEBUG_STDOUT
  groupcount = clpPathGroupVector_GetCount( &clpoffset->grouplist );
  for( groupindex = 0 ; groupindex < groupcount ; groupindex++ )
    clpDebugShape( &(clpPathGroupVector_Access( &clpoffset->grouplist, groupindex )->inputshape), "Shape" );
#endif

#if CLP_OFFSET_DEBUG_IMAGES
  for( groupindex = 0 ; groupindex < groupcount ; groupindex++ )
    clpDebugDrawShape( &debugimage, &(clpPathGroupVector_Access( &clpoffset->grouplist, groupindex )->inputshape), (float[]){ 128.0f, 128.0f, 0.0f, 255.0f } );
  snprintf( debugimagepath, 256, "debug-polyoffset/%05db.png", debugimageindex );
  clpDebugSaveImage( &debugimage, debugimagepath );
#endif

  if( clpoffset->mergegroupsflag && ( groupcount > 0 ) )
  {
    reversedflag = clpPathGroupVector_Access( &clpoffset->grouplist, 0 )->reversedflag;

    /* Clean up self-intersections */
    clp = clpCreateClipper( 0 );

    clpSetPreserveCollinear( clp, 0 );
    /* The solution should retain the orientation of the input */
    clpSetReverseSolution( clp, ( clpoffset->reversesolutionflag != reversedflag ) );

    workshape.pathcount = clpPathVector_GetCount( &clpoffset->pathlist );
    workshape.pathlist = clpPathVector_Steal( &clpoffset->pathlist );

#if CLP_OFFSET_DEBUG_STDOUT
    clpDebugShape( &workshape, "Offset" );
#endif

#if CLP_OFFSET_DEBUG_IMAGES
    clpDebugDrawShape( &debugimage, &workshape, (float[]){ 0.0f, 0.0f, 255.0f, 255.0f } );
#endif

    clpAddShape( clp, &workshape, CLP_PATHTYPE_SUBJECT, 0 );

    if( reversedflag )
      clpClipperExecute( clp, CLP_CLIPTYPE_UNION, CLP_FILLRULE_NEGATIVE, &clpoffset->finalshape );
    else
      clpClipperExecute( clp, CLP_CLIPTYPE_UNION, CLP_FILLRULE_POSITIVE, &clpoffset->finalshape );

#if CLP_OFFSET_DEBUG_STDOUT
    clpDebugShape( &clpoffset->finalshape, "Final" );
#endif

#if CLP_OFFSET_DEBUG_IMAGES
    clpDebugDrawShape( &debugimage, &clpoffset->finalshape, (float[]){ 128.0f, 255.0f, 255.0f, 255.0f } );
#endif

    clpFreeShape( &workshape );
    clpDestroyClipper( clp );
  }
  else
  {
    clpoffset->finalshape.pathcount = clpPathVector_GetCount( &clpoffset->pathlist );
    clpoffset->finalshape.pathlist = clpPathVector_Steal( &clpoffset->pathlist );
  }

#if CLP_OFFSET_DEBUG_IMAGES
  snprintf( debugimagepath, 256, "debug-polyoffset/%05dc.png", debugimageindex );
  clpDebugSaveImage( &debugimage, debugimagepath );
  clpDebugFreeImage( &debugimage );
  debugimageindex++;
#endif

  return clpoffset->finalshape;
}



clpOffset *clpOffsetCreate()
{
  clpOffset *clpoffset;
  clpoffset = (clpOffset *)malloc( sizeof(clpOffset) );
  memset( (void *)clpoffset, 0, sizeof(clpOffset) );
  /* 16 steps for full circle */
  clpoffset->roundradiansteps = 16.0 / ( 2 * M_PI );
  clpoffset->miterlimit = 2.0;
  clpoffset->preservecollinearflag = 0;
  clpoffset->reversesolutionflag = 0;
  clpoffset->jointype = CLP_JOINTYPE_SQUARE;
  clpoffset->mergegroupsflag = 1;
  clpPathD_Init( &clpoffset->norms, 0 );
  return clpoffset;
}


void clpOffsetDestroy( clpOffset *clpoffset )
{
  size_t groupindex, groupcount;
  clpPathGroup *pathgroup;

  groupcount = clpPathGroupVector_GetCount( &clpoffset->grouplist );
  for( groupindex = 0 ; groupindex < groupcount ; groupindex++ )
  {
    pathgroup = clpPathGroupVector_Access( &clpoffset->grouplist, groupindex );
    if( pathgroup->inputpathlistallocflag )
      free( (pathgroup->inputshape).pathlist );
  }
  clpPathVector_Clear( &clpoffset->pathlist );
  clpPathGroupVector_Clear( &clpoffset->grouplist );
  clpPathD_Clear( &clpoffset->norms );
  free( clpoffset );
  return;
}