static VALUE cmethod_create(VALUE module, VALUE factory, VALUE exterior, VALUE interior_array)
{
  Check_Type(interior_array, T_ARRAY);
  RGeo_FactoryData* factory_data = RGEO_FACTORY_DATA_PTR(factory);
  VALUE linear_ring_type = factory_data->globals->feature_linear_ring;
  GEOSGeometry* exterior_geom = rgeo_convert_to_detached_geos_geometry(exterior, factory, linear_ring_type, NULL);
  if (exterior_geom) {
    GEOSContextHandle_t context = factory_data->geos_context;
    unsigned int len = (unsigned int)RARRAY_LEN(interior_array);
    GEOSGeometry** interior_geoms = ALLOC_N(GEOSGeometry*, len == 0 ? 1 : len);
    if (interior_geoms) {
      unsigned int actual_len = 0;
      unsigned int i;
      for (i=0; i<len; ++i) {
        GEOSGeometry* interior_geom = rgeo_convert_to_detached_geos_geometry(rb_ary_entry(interior_array, i), factory, linear_ring_type, NULL);
        if (interior_geom) {
          interior_geoms[actual_len++] = interior_geom;
        }
      }
      if (len == actual_len) {
        GEOSGeometry* polygon = GEOSGeom_createPolygon_r(context, exterior_geom, interior_geoms, actual_len);
        if (polygon) {
          free(interior_geoms);
          return rgeo_wrap_geos_geometry(factory, polygon, factory_data->globals->geos_polygon);
        }
      }
      for (i=0; i<actual_len; ++i) {
        GEOSGeom_destroy_r(context, interior_geoms[i]);
      }
      free(interior_geoms);
    }
    GEOSGeom_destroy_r(context, exterior_geom);
  }
  return Qnil;
}

// Return closest point to given distance within geometry.
// 'spgeom' must be a LineString
SEXP rgeos_interpolate(SEXP env, SEXP spgeom, SEXP d, SEXP normalized) {

    GEOSContextHandle_t GEOShandle = getContextHandle(env);
    GEOSGeom geom = rgeos_convert_R2geos(env, spgeom);

    GEOSGeom res_geos;
    double dist;

    int nlines = length(GET_SLOT(spgeom, install("lines")));
    if (nlines < 1) {
        error("rgeos_project: invalid number of lines");
    }

    int n = LENGTH(d);
    if (n < 1) {
        error("rgeos_interpolate: invalid number of requested points");
    }

    int pc = 0;
    SEXP crd;
    PROTECT(crd = NEW_NUMERIC(n*2)); pc++;

    double x;
    double y;
    SEXP ans;

    // select interpolation function (normalized/unnormalized)
    GEOSGeometry GEOS_DLL *(*interp_fun)(GEOSContextHandle_t,
                                         const GEOSGeometry*,
                                         double);

    if (LOGICAL_POINTER(normalized)[0]) {
        interp_fun = &GEOSInterpolateNormalized_r;
    } else {
        interp_fun = &GEOSInterpolate_r;
    }

    // interpolate points and store result in coord matrix
    for (int i = 0; i < n; i++) {

        dist = NUMERIC_POINTER(d)[i];

        res_geos = (*interp_fun)(GEOShandle, geom, dist);

        rgeos_Pt2xy(env, res_geos, &x, &y);

        NUMERIC_POINTER(crd)[i] = x;
        NUMERIC_POINTER(crd)[n+i] = y;
    }

    GEOSGeom_destroy_r(GEOShandle, geom);
    GEOSGeom_destroy_r(GEOShandle, res_geos);

    // return coordinates as matrix
    PROTECT(ans = rgeos_formatcrdMat(crd, n)); pc++;

    UNPROTECT(pc);
    return(ans);
}

void QgsZonalStatistics::statisticsFromMiddlePointTest( void* band, const QgsGeometry& poly, int pixelOffsetX,
    int pixelOffsetY, int nCellsX, int nCellsY, double cellSizeX, double cellSizeY, const QgsRectangle& rasterBBox, FeatureStats &stats )
{
  double cellCenterX, cellCenterY;

  float* scanLine = ( float * ) CPLMalloc( sizeof( float ) * nCellsX );
  cellCenterY = rasterBBox.yMaximum() - pixelOffsetY * cellSizeY - cellSizeY / 2;
  stats.reset();

  GEOSGeometry* polyGeos = poly.exportToGeos();
  if ( !polyGeos )
  {
    return;
  }

  GEOSContextHandle_t geosctxt = QgsGeometry::getGEOSHandler();
  const GEOSPreparedGeometry* polyGeosPrepared = GEOSPrepare_r( geosctxt, polyGeos );
  if ( !polyGeosPrepared )
  {
    GEOSGeom_destroy_r( geosctxt, polyGeos );
    return;
  }

  GEOSCoordSequence* cellCenterCoords = nullptr;
  GEOSGeometry* currentCellCenter = nullptr;

  for ( int i = 0; i < nCellsY; ++i )
  {
    if ( GDALRasterIO( band, GF_Read, pixelOffsetX, pixelOffsetY + i, nCellsX, 1, scanLine, nCellsX, 1, GDT_Float32, 0, 0 )
         != CPLE_None )
    {
      continue;
    }
    cellCenterX = rasterBBox.xMinimum() + pixelOffsetX * cellSizeX + cellSizeX / 2;
    for ( int j = 0; j < nCellsX; ++j )
    {
      if ( validPixel( scanLine[j] ) )
      {
        GEOSGeom_destroy_r( geosctxt, currentCellCenter );
        cellCenterCoords = GEOSCoordSeq_create_r( geosctxt, 1, 2 );
        GEOSCoordSeq_setX_r( geosctxt, cellCenterCoords, 0, cellCenterX );
        GEOSCoordSeq_setY_r( geosctxt, cellCenterCoords, 0, cellCenterY );
        currentCellCenter = GEOSGeom_createPoint_r( geosctxt, cellCenterCoords );
        if ( GEOSPreparedContains_r( geosctxt, polyGeosPrepared, currentCellCenter ) )
        {
          stats.addValue( scanLine[j] );
        }
      }
      cellCenterX += cellSizeX;
    }
    cellCenterY -= cellSizeY;
  }
  GEOSGeom_destroy_r( geosctxt, currentCellCenter );
  CPLFree( scanLine );
  GEOSPreparedGeom_destroy_r( geosctxt, polyGeosPrepared );
  GEOSGeom_destroy_r( geosctxt, polyGeos );
}

void QgsZonalStatistics::statisticsFromMiddlePointTest( const QgsGeometry &poly, int pixelOffsetX,
    int pixelOffsetY, int nCellsX, int nCellsY, double cellSizeX, double cellSizeY, const QgsRectangle &rasterBBox, FeatureStats &stats )
{
  double cellCenterX, cellCenterY;

  cellCenterY = rasterBBox.yMaximum() - pixelOffsetY * cellSizeY - cellSizeY / 2;
  stats.reset();

  GEOSGeometry *polyGeos = poly.exportToGeos();
  if ( !polyGeos )
  {
    return;
  }

  GEOSContextHandle_t geosctxt = QgsGeometry::getGEOSHandler();
  const GEOSPreparedGeometry *polyGeosPrepared = GEOSPrepare_r( geosctxt, polyGeos );
  if ( !polyGeosPrepared )
  {
    GEOSGeom_destroy_r( geosctxt, polyGeos );
    return;
  }

  GEOSCoordSequence *cellCenterCoords = nullptr;
  GEOSGeometry *currentCellCenter = nullptr;

  QgsRectangle featureBBox = poly.boundingBox().intersect( &rasterBBox );
  QgsRectangle intersectBBox = rasterBBox.intersect( &featureBBox );

  QgsRasterBlock *block = mRasterProvider->block( mRasterBand, intersectBBox, nCellsX, nCellsY );
  for ( int i = 0; i < nCellsY; ++i )
  {
    cellCenterX = rasterBBox.xMinimum() + pixelOffsetX * cellSizeX + cellSizeX / 2;
    for ( int j = 0; j < nCellsX; ++j )
    {
      if ( validPixel( block->value( i, j ) ) )
      {
        GEOSGeom_destroy_r( geosctxt, currentCellCenter );
        cellCenterCoords = GEOSCoordSeq_create_r( geosctxt, 1, 2 );
        GEOSCoordSeq_setX_r( geosctxt, cellCenterCoords, 0, cellCenterX );
        GEOSCoordSeq_setY_r( geosctxt, cellCenterCoords, 0, cellCenterY );
        currentCellCenter = GEOSGeom_createPoint_r( geosctxt, cellCenterCoords );
        if ( GEOSPreparedContains_r( geosctxt, polyGeosPrepared, currentCellCenter ) )
        {
          stats.addValue( block->value( i, j ) );
        }
      }
      cellCenterX += cellSizeX;
    }
    cellCenterY -= cellSizeY;
  }

  GEOSGeom_destroy_r( geosctxt, currentCellCenter );
  GEOSPreparedGeom_destroy_r( geosctxt, polyGeosPrepared );
  GEOSGeom_destroy_r( geosctxt, polyGeos );
  delete block;
}

QgsLabelFeature::~QgsLabelFeature()
{
    if ( mGeometry )
        GEOSGeom_destroy_r( QgsGeometry::getGEOSHandler(), mGeometry );

    if ( mObstacleGeometry )
        GEOSGeom_destroy_r( QgsGeometry::getGEOSHandler(), mObstacleGeometry );

    delete mInfo;
}

OGRMultiPolygon* OGRILI1Layer::Polygonize( OGRGeometryCollection* poLines, bool fix_crossing_lines )
{
    OGRMultiPolygon *poPolygon = new OGRMultiPolygon();

    if (poLines->getNumGeometries() == 0) return poPolygon;

#if defined(HAVE_GEOS)
    GEOSGeom *ahInGeoms = NULL;
    int       i = 0;
    OGRGeometryCollection *poNoncrossingLines = poLines;
    GEOSGeom hResultGeom = NULL;
    OGRGeometry *poMP = NULL;

    if (fix_crossing_lines && poLines->getNumGeometries() > 0)
    {
        CPLDebug( "OGR_ILI", "Fixing crossing lines");
        //A union of the geometry collection with one line fixes invalid geometries
        poNoncrossingLines = (OGRGeometryCollection*)poLines->Union(poLines->getGeometryRef(0));
        CPLDebug( "OGR_ILI", "Fixed lines: %d", poNoncrossingLines->getNumGeometries()-poLines->getNumGeometries());
    }
    
    GEOSContextHandle_t hGEOSCtxt = OGRGeometry::createGEOSContext();

    ahInGeoms = (GEOSGeom *) CPLCalloc(sizeof(void*),poNoncrossingLines->getNumGeometries());
    for( i = 0; i < poNoncrossingLines->getNumGeometries(); i++ )
          ahInGeoms[i] = poNoncrossingLines->getGeometryRef(i)->exportToGEOS(hGEOSCtxt);

    hResultGeom = GEOSPolygonize_r( hGEOSCtxt,
                                    ahInGeoms,
                                   poNoncrossingLines->getNumGeometries() );

    for( i = 0; i < poNoncrossingLines->getNumGeometries(); i++ )
        GEOSGeom_destroy_r( hGEOSCtxt, ahInGeoms[i] );
    CPLFree( ahInGeoms );
    if (poNoncrossingLines != poLines) delete poNoncrossingLines;

    if( hResultGeom == NULL )
    {
        OGRGeometry::freeGEOSContext( hGEOSCtxt );
        return NULL;
    }

    poMP = OGRGeometryFactory::createFromGEOS( hGEOSCtxt, hResultGeom );

    GEOSGeom_destroy_r( hGEOSCtxt, hResultGeom );
    OGRGeometry::freeGEOSContext( hGEOSCtxt );

    return (OGRMultiPolygon *) poMP;

#endif

    return poPolygon;
}

// Return distance of points 'spppoints' projected on 'spgeom' from origin
// of 'spgeom'. Geometry 'spgeom' must be a lineal geometry
SEXP rgeos_project(SEXP env, SEXP spgeom, SEXP sppoint, SEXP normalized) {

    GEOSContextHandle_t GEOShandle = getContextHandle(env);
    GEOSGeom geom = rgeos_convert_R2geos(env, spgeom);

    SEXP crds = GET_SLOT(sppoint, install("coords"));
    SEXP dim = getAttrib(crds, install("dim"));

    int nlines = length(GET_SLOT(spgeom, install("lines")));
    if (nlines < 1) {
        error("rgeos_project: invalid number of lines");
    }

    int n = INTEGER_POINTER(dim)[0];
    if (n < 1) {
        error("rgeos_project: invalid number of points");
    }

    int pc = 0;
    SEXP ans;
    PROTECT(ans = NEW_NUMERIC(n)); pc++;

    GEOSGeom p;

    // select projection function (normalized/unnormalized)
    double GEOS_DLL (*proj_fun)(GEOSContextHandle_t,
                                const GEOSGeometry*,
                                const GEOSGeometry*);

    if (LOGICAL_POINTER(normalized)[0]) {
        proj_fun = &GEOSProjectNormalized_r;
    } else {
        proj_fun = &GEOSProject_r;
    }

    // project points to line geometry
    for (int i = 0; i < n; i++) {

        p = rgeos_xy2Pt(env,
                        NUMERIC_POINTER(crds)[i],
                        NUMERIC_POINTER(crds)[i+n]);

        NUMERIC_POINTER(ans)[i] = (*proj_fun)(GEOShandle, geom, p);
    }

    GEOSGeom_destroy_r(GEOShandle, geom);
    GEOSGeom_destroy_r(GEOShandle, p);

    UNPROTECT(pc);

    return(ans);
}

SEXP rgeos_node(SEXP env, SEXP obj) {

    SEXP ans, id;
    int pc=0;
    
    GEOSContextHandle_t GEOShandle = getContextHandle(env);
    SEXP p4s = GET_SLOT(obj, install("proj4string"));
    GEOSGeom geom = rgeos_convert_R2geos(env, obj);
//    int type = GEOSGeomTypeId_r(GEOShandle, geom);
//Rprintf("type: %d, %s\n", type, GEOSGeomType_r(GEOShandle, geom));
    
    GEOSGeom res = GEOSNode_r(GEOShandle, geom);
    
//    type = GEOSGeomTypeId_r(GEOShandle, res);

    int ng = GEOSGetNumGeometries_r(GEOShandle, res);

//Rprintf("ng: %d, type: %d, %s\n", ng, type, GEOSGeomType_r(GEOShandle, res));

    char buf[BUFSIZ];

    PROTECT(id = NEW_CHARACTER(ng)); pc++;
    for (int i=0; i<ng; i++) {
        sprintf(buf, "%d", i);
        SET_STRING_ELT(id, i, COPY_TO_USER_STRING(buf));
    }

    GEOSGeom_destroy_r(GEOShandle, geom);

    ans = rgeos_convert_geos2R(env, res, p4s, id); 

    UNPROTECT(pc);
    return(ans);

}

static VALUE method_geometry_envelope(VALUE self)
{
  VALUE result;
  RGeo_GeometryData* self_data;
  const GEOSGeometry* self_geom;
  GEOSContextHandle_t geos_context;
  GEOSGeometry* envelope;

  result = Qnil;
  self_data = RGEO_GEOMETRY_DATA_PTR(self);
  self_geom = self_data->geom;
  if (self_geom) {
    geos_context = self_data->geos_context;
    envelope = GEOSEnvelope_r(geos_context, self_geom);
    // GEOS returns an "empty" point for an empty collection's envelope.
    // We don't allow that type, so we replace it with an empty collection.
    if (!envelope ||
        GEOSGeomTypeId_r(geos_context, envelope) == GEOS_POINT &&
        GEOSGetNumCoordinates_r(geos_context, envelope) == 0) {
      if (envelope) {
        GEOSGeom_destroy_r(geos_context, envelope);
      }
      envelope = GEOSGeom_createCollection_r(geos_context, GEOS_GEOMETRYCOLLECTION, NULL, 0);
    }
    result = rgeo_wrap_geos_geometry(self_data->factory, envelope, Qnil);
  }
  return result;
}

static VALUE method_geometry_initialize_copy(VALUE self, VALUE orig)
{
  // Clear out any existing value
  RGeo_GeometryData* self_data = RGEO_GEOMETRY_DATA_PTR(self);
  GEOSGeometry* self_geom = self_data->geom;
  if (self_geom) {
    GEOSGeom_destroy_r(self_data->geos_context, self_geom);
    self_data->geom = NULL;
    self_data->geos_context = NULL;
    self_data->factory = Qnil;
    self_data->klasses = Qnil;
  }
  
  // Copy value from orig
  const GEOSGeometry* geom = rgeo_get_geos_geometry_safe(orig);
  if (geom) {
    RGeo_GeometryData* orig_data = RGEO_GEOMETRY_DATA_PTR(orig);
    GEOSContextHandle_t orig_context = orig_data->geos_context;
    GEOSGeometry* clone_geom = GEOSGeom_clone_r(orig_context, geom);
    if (clone_geom) {
      GEOSSetSRID_r(orig_context, clone_geom, GEOSGetSRID_r(orig_context, geom));
      self_data->geom = clone_geom;
      self_data->geos_context = orig_context;
      self_data->factory = orig_data->factory;
      self_data->klasses = orig_data->klasses;
    }
  }
  return self;
}

BOX3D Polygon::bounds() const
{
    uint32_t numInputDims;
    BOX3D output;

    GEOSGeometry* boundary = GEOSGeom_clone_r(m_ctx, m_geom);

    // Smash out multi*
    if (GEOSGeomTypeId_r(m_ctx, m_geom) > 3)
        boundary = GEOSEnvelope_r(m_ctx, m_geom);

    GEOSGeometry const* ring = GEOSGetExteriorRing_r(m_ctx, boundary);
    GEOSCoordSequence const* coords = GEOSGeom_getCoordSeq_r(m_ctx, ring);

    GEOSCoordSeq_getDimensions_r(m_ctx, coords, &numInputDims);

    uint32_t count(0);
    GEOSCoordSeq_getSize_r(m_ctx, coords, &count);

    double x(0.0);
    double y(0.0);
    double z(0.0);
    for (unsigned i = 0; i < count; ++i)
    {
        GEOSCoordSeq_getOrdinate_r(m_ctx, coords, i, 0, &x);
        GEOSCoordSeq_getOrdinate_r(m_ctx, coords, i, 1, &y);
        if (numInputDims > 2)
            GEOSCoordSeq_getOrdinate_r(m_ctx, coords, i, 2, &z);
        output.grow(x, y, z);
    }
    GEOSGeom_destroy_r(m_ctx, boundary);

    return output;
}

SEXP rgeos_miscfunc(SEXP env, SEXP obj, SEXP byid, p_miscfunc miscfunc) {

    SEXP ans;
    
    GEOSContextHandle_t GEOShandle = getContextHandle(env);

    GEOSGeom geom = rgeos_convert_R2geos(env, obj);
    int type = GEOSGeomTypeId_r(GEOShandle, geom);
    
    int n = (LOGICAL_POINTER(byid)[0] && type == GEOS_GEOMETRYCOLLECTION) ? 
                GEOSGetNumGeometries_r(GEOShandle, geom) : 1;
    
    int pc=0;
    PROTECT(ans = NEW_NUMERIC(n)); pc++;

    GEOSGeom curgeom = geom;
    for(int i=0; i<n; i++) {
        if ( n > 1) {
            curgeom = (GEOSGeom) GEOSGetGeometryN_r(GEOShandle, geom, i);
            if (curgeom == NULL) error("rgeos_miscfunc: unable to get subgeometries");
        }
        
        double val;
        if (!miscfunc(GEOShandle, curgeom, &val))
            error("rgeos_miscfunc: unable to calculate");
            
        NUMERIC_POINTER(ans)[i] = val;
    }

    GEOSGeom_destroy_r(GEOShandle, geom);

    UNPROTECT(pc);
    return(ans);
}

// Fully node given linework
static GEOSGeometry *LWGEOM_GEOS_nodeLines( const GEOSGeometry *lines )
{
  GEOSContextHandle_t handle = QgsGeos::getGEOSHandler();

  // Union with first geometry point, obtaining full noding
  // and dissolving of duplicated repeated points
  //
  // TODO: substitute this with UnaryUnion?

  GEOSGeometry *point = LWGEOM_GEOS_getPointN( lines, 0 );
  if ( ! point )
    return nullptr;

  GEOSGeometry *noded = nullptr;
  try
  {
    noded = GEOSUnion_r( handle, lines, point );
  }
  catch ( GEOSException & )
  {
    // no need to do anything here - we'll return nullptr anyway
  }
  GEOSGeom_destroy_r( handle, point );
  return noded;
}

bool Polygon::covers(PointRef& ref) const
{
    GEOSCoordSequence* coords = GEOSCoordSeq_create_r(m_ctx, 1, 3);
    if (!coords)
        throw pdal_error("Unable to allocate coordinate sequence");

    const double x = ref.getFieldAs<double>(Dimension::Id::X);
    const double y = ref.getFieldAs<double>(Dimension::Id::Y);
    const double z = ref.getFieldAs<double>(Dimension::Id::Z);

    if (!GEOSCoordSeq_setX_r(m_ctx, coords, 0, x))
        throw pdal_error("unable to set x for coordinate sequence");
    if (!GEOSCoordSeq_setY_r(m_ctx, coords, 0, y))
        throw pdal_error("unable to set y for coordinate sequence");
    if (!GEOSCoordSeq_setZ_r(m_ctx, coords, 0, z))
        throw pdal_error("unable to set z for coordinate sequence");
    GEOSGeometry* p = GEOSGeom_createPoint_r(m_ctx, coords);
    if (!p)
        throw pdal_error("unable to allocate candidate test point");

    bool covers = (bool)(GEOSPreparedCovers_r(m_ctx, m_prepGeom, p));
    GEOSGeom_destroy_r(m_ctx, p);

    return covers;
}

  int LabelPosition::polygonIntersectionCost( PointSet *polygon ) const
  {
    if ( !mGeos )
      createGeosGeom();

    if ( !polygon->mGeos )
      polygon->createGeosGeom();

    GEOSContextHandle_t geosctxt = geosContext();

    int cost = 0;
    //check the label center. if covered by polygon, initial cost of 4
    if ( polygon->containsPoint(( x[0] + x[2] ) / 2.0, ( y[0] + y[2] ) / 2.0 ) )
      cost += 4;

    try
    {
      //calculate proportion of label candidate which is covered by polygon
      GEOSGeometry* intersectionGeom = GEOSIntersection_r( geosctxt, mGeos, polygon->mGeos );
      if ( !intersectionGeom )
        return cost;

      double positionArea = 0;
      if ( GEOSArea_r( geosctxt, mGeos, &positionArea ) != 1 )
      {
        GEOSGeom_destroy_r( geosctxt, intersectionGeom );
        return cost;
      }

      double intersectionArea = 0;
      if ( GEOSArea_r( geosctxt, intersectionGeom, &intersectionArea ) != 1 )
      {
        intersectionArea = 0;
      }

      GEOSGeom_destroy_r( geosctxt, intersectionGeom );

      double portionCovered = intersectionArea / positionArea;
      cost += ceil( portionCovered * 8.0 ); //cost of 8 if totally covered
      return cost;
    }
    catch ( GEOSException &e )
    {
      QgsMessageLog::logMessage( QObject::tr( "Exception: %1" ).arg( e.what() ), QObject::tr( "GEOS" ) );
      return cost;
    }
  }

void Crop::crop(PointBuffer& input, PointBuffer& output)
{
    bool logOutput = (log()->getLevel() > LogLevel::Debug4);
    if (logOutput)
        log()->floatPrecision(8);

    for (PointId idx = 0; idx < input.size(); ++idx)
    {
        double x = input.getFieldAs<double>(Dimension::Id::X, idx);
        double y = input.getFieldAs<double>(Dimension::Id::Y, idx);
        double z = input.getFieldAs<double>(Dimension::Id::Z, idx);

        if (logOutput)
        {
            log()->floatPrecision(10);
            log()->get(LogLevel::Debug5) << "input: " << x << " y: " << y <<
                " z: " << z << std::endl;
        }

        if (m_poly.empty())
        {
            // We don't have a polygon, just a bounds. Filter on that
            // by itself.
            if (!m_cropOutside && m_bounds.contains(x, y, z))
                output.appendPoint(input, idx);
        }
#ifdef PDAL_HAVE_GEOS
        else
        {
            int ret(0);

            // precise filtering based on the geometry
            GEOSCoordSequence* coords =
                GEOSCoordSeq_create_r(m_geosEnvironment, 1, 3);
            if (!coords)
                throw pdal_error("unable to allocate coordinate sequence");
            ret = GEOSCoordSeq_setX_r(m_geosEnvironment, coords, 0, x);
            if (!ret)
                throw pdal_error("unable to set x for coordinate sequence");
            ret = GEOSCoordSeq_setY_r(m_geosEnvironment, coords, 0, y);
            if (!ret)
                throw pdal_error("unable to set y for coordinate sequence");
            ret = GEOSCoordSeq_setZ_r(m_geosEnvironment, coords, 0, z);
            if (!ret)
                throw pdal_error("unable to set z for coordinate sequence");

            GEOSGeometry* p = GEOSGeom_createPoint_r(m_geosEnvironment, coords);
            if (!p)
                throw pdal_error("unable to allocate candidate test point");

            if (static_cast<bool>(GEOSPreparedContains_r(m_geosEnvironment,
                m_geosPreparedGeometry, p)) != m_cropOutside)
                output.appendPoint(input, idx);
            GEOSGeom_destroy_r(m_geosEnvironment, p);
        }
#endif
    }
}

Polygon::~Polygon()
{
    if (m_geom)
        GEOSGeom_destroy_r(m_ctx, m_geom);
    if (m_prepGeom)
        GEOSPreparedGeom_destroy_r(m_ctx, m_prepGeom);
    m_geom = 0;
    m_prepGeom = 0;
}

static bool _canExportToGeos( const QgsGeometry &geom )
{
  GEOSGeometry *geosGeom = geom.exportToGeos();
  if ( geosGeom )
  {
    GEOSGeom_destroy_r( QgsGeometry::getGEOSHandler(), geosGeom );
    return true;
  }
  return false;
}

void Crop::done(PointContext ctx)
{
#ifdef PDAL_HAVE_GEOS
    if (m_geosPreparedGeometry)
        GEOSPreparedGeom_destroy_r(m_geosEnvironment, m_geosPreparedGeometry);

    if (m_geosGeometry)
        GEOSGeom_destroy_r(m_geosEnvironment, m_geosGeometry);

    if (m_geosEnvironment)
        finishGEOS_r(m_geosEnvironment);
#endif
}

static void free_geos_geom(void* data) {
  if (data == NULL) {
    return;
  }

  geos_geometry_t* geom = (geos_geometry_t*)data;
  GEOSGeom_destroy_r(geom->context, geom->geometry);

  geom->context = NULL;
  geom->geometry = NULL;

  sqlite3_free(data);
}