/******************************************************************************
 * $Id$
 *
 * Project:  MapServer
 * Purpose:  MapServer-GEOS integration.
 * Author:   Steve Lime and the MapServer team.
 *
 ******************************************************************************
 * Copyright (c) 1996-2005 Regents of the University of Minnesota.
 *
 * Permission is hereby granted, free of charge, to any person obtaining a
 * copy of this software and associated documentation files (the "Software"),
 * to deal in the Software without restriction, including without limitation
 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
 * and/or sell copies of the Software, and to permit persons to whom the
 * Software is furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in
 * all copies of this Software or works derived from this Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
 * DEALINGS IN THE SOFTWARE.
 *****************************************************************************/

#include "mapserver.h"
#include "mapgraph.h"

#ifdef USE_GEOS

// To avoid accidental use of non reentrant GEOS API.
// (check only effective in GEOS >= 3.5)
#define GEOS_USE_ONLY_R_API

#include <geos_c.h>

/*
** Error handling...
*/
static void msGEOSError(const char *format, ...) {
  va_list args;

  va_start(args, format);
  msSetError(MS_GEOSERR, format, "msGEOSError()",
             args); /* just pass along to MapServer error handling */
  va_end(args);

  return;
}

static void msGEOSNotice(const char *fmt, ...) {
  (void)fmt; /* do nothing with notices at this point */
}

#ifndef USE_THREAD
static GEOSContextHandle_t geos_handle;

static inline GEOSContextHandle_t msGetGeosContextHandle() {
  return geos_handle;
}

#else
#include "mapthread.h"
typedef struct geos_thread_info {
  struct geos_thread_info *next;
  void *thread_id;
  GEOSContextHandle_t geos_handle;
} geos_thread_info_t;

static geos_thread_info_t *geos_list = NULL;

static GEOSContextHandle_t msGetGeosContextHandle() {
  geos_thread_info_t *link;
  GEOSContextHandle_t ret_obj;
  void *thread_id;

  msAcquireLock(TLOCK_GEOS);

  thread_id = msGetThreadId();

  /* find link for this thread */

  for (link = geos_list;
       link != NULL && link->thread_id != thread_id && link->next != NULL &&
       link->next->thread_id != thread_id;
       link = link->next) {
  }

  /* If the target thread link is already at the head of the list were ok */
  if (geos_list != NULL && geos_list->thread_id == thread_id) {
  }

  /* We don't have one ... initialize one. */
  else if (link == NULL || link->next == NULL) {
    geos_thread_info_t *new_link;
    new_link = (geos_thread_info_t *)malloc(sizeof(geos_thread_info_t));
    new_link->next = geos_list;
    new_link->thread_id = thread_id;
    new_link->geos_handle = initGEOS_r(msGEOSNotice, msGEOSError);

    geos_list = new_link;
  }

  /* If the link is not already at the head of the list, promote it */
  else {
    geos_thread_info_t *target = link->next;

    link->next = link->next->next;
    target->next = geos_list;
    geos_list = target;
  }

  ret_obj = geos_list->geos_handle;

  msReleaseLock(TLOCK_GEOS);

  return ret_obj;
}
#endif

/*
** Setup/Cleanup wrappers
*/
void msGEOSSetup() {
#ifndef USE_THREAD
  geos_handle = initGEOS_r(msGEOSNotice, msGEOSError);
#endif
}

void msGEOSCleanup() {
#ifndef USE_THREAD
  finishGEOS_r(geos_handle);
  geos_handle = NULL;
#else
  geos_thread_info_t *link;
  msAcquireLock(TLOCK_GEOS);
  for (link = geos_list; link != NULL;) {
    geos_thread_info_t *cur = link;
    link = link->next;
    finishGEOS_r(cur->geos_handle);
    free(cur);
  }
  geos_list = NULL;
  msReleaseLock(TLOCK_GEOS);
#endif
}

/*
** Translation functions
*/
static GEOSGeom msGEOSShape2Geometry_point(pointObj *point) {
  GEOSCoordSeq coords;
  GEOSGeom g;
  GEOSContextHandle_t handle = msGetGeosContextHandle();

  if (!point)
    return NULL;

  coords = GEOSCoordSeq_create_r(handle, 1, 2); /* todo handle z's */
  if (!coords)
    return NULL;

  GEOSCoordSeq_setX_r(handle, coords, 0, point->x);
  GEOSCoordSeq_setY_r(handle, coords, 0, point->y);
  /* GEOSCoordSeq_setY(coords, 0, point->z); */

  g = GEOSGeom_createPoint_r(handle,
                             coords); /* g owns the coordinate in coords */

  return g;
}

static GEOSGeom msGEOSShape2Geometry_multipoint(lineObj *multipoint) {
  int i;
  GEOSGeom g;
  GEOSGeom *points;
  GEOSContextHandle_t handle = msGetGeosContextHandle();

  if (!multipoint)
    return NULL;

  points = malloc(multipoint->numpoints * sizeof(GEOSGeom));
  if (!points)
    return NULL;

  for (i = 0; i < multipoint->numpoints; i++)
    points[i] = msGEOSShape2Geometry_point(&(multipoint->point[i]));

  g = GEOSGeom_createCollection_r(handle, GEOS_MULTIPOINT, points,
                                  multipoint->numpoints);

  free(points);

  return g;
}

static GEOSGeom msGEOSShape2Geometry_line(lineObj *line) {
  int i;
  GEOSGeom g;
  GEOSCoordSeq coords;
  GEOSContextHandle_t handle = msGetGeosContextHandle();

  if (!line)
    return NULL;

  coords =
      GEOSCoordSeq_create_r(handle, line->numpoints, 2); /* todo handle z's */
  if (!coords)
    return NULL;

  for (i = 0; i < line->numpoints; i++) {
    GEOSCoordSeq_setX_r(handle, coords, i, line->point[i].x);
    GEOSCoordSeq_setY_r(handle, coords, i, line->point[i].y);
    /* GEOSCoordSeq_setZ(coords, i, line->point[i].z); */
  }

  g = GEOSGeom_createLineString_r(
      handle, coords); /* g owns the coordinates in coords */

  return g;
}

static GEOSGeom msGEOSShape2Geometry_multiline(shapeObj *multiline) {
  int i;
  GEOSGeom g;
  GEOSGeom *lines;
  GEOSContextHandle_t handle = msGetGeosContextHandle();

  if (!multiline)
    return NULL;

  lines = malloc(multiline->numlines * sizeof(GEOSGeom));
  if (!lines)
    return NULL;

  for (i = 0; i < multiline->numlines; i++)
    lines[i] = msGEOSShape2Geometry_line(&(multiline->line[i]));

  g = GEOSGeom_createCollection_r(handle, GEOS_MULTILINESTRING, lines,
                                  multiline->numlines);

  free(lines);

  return g;
}

static GEOSGeom msGEOSShape2Geometry_simplepolygon(shapeObj *shape, int r,
                                                   int *outerList) {
  int i, j, k;
  GEOSCoordSeq coords;
  GEOSGeom g;
  GEOSGeom outerRing;
  GEOSGeom *innerRings = NULL;
  int numInnerRings = 0, *innerList;
  GEOSContextHandle_t handle = msGetGeosContextHandle();

  if (!shape || !outerList)
    return NULL;

  /* build the outer shell */
  coords = GEOSCoordSeq_create_r(handle, shape->line[r].numpoints,
                                 2); /* todo handle z's */
  if (!coords)
    return NULL;

  for (i = 0; i < shape->line[r].numpoints; i++) {
    GEOSCoordSeq_setX_r(handle, coords, i, shape->line[r].point[i].x);
    GEOSCoordSeq_setY_r(handle, coords, i, shape->line[r].point[i].y);
    /* GEOSCoordSeq_setZ(coords, i, shape->line[r].point[i].z); */
  }

  outerRing = GEOSGeom_createLinearRing_r(
      handle, coords); /* outerRing owns the coordinates in coords */

  /* build the holes */
  innerList = msGetInnerList(shape, r, outerList);
  for (j = 0; j < shape->numlines; j++)
    if (innerList[j] == MS_TRUE)
      numInnerRings++;

  if (numInnerRings > 0) {
    k = 0; /* inner ring counter */

    innerRings = msSmallMalloc(numInnerRings * sizeof(GEOSGeom));

    for (j = 0; j < shape->numlines; j++) {
      if (innerList[j] == MS_FALSE)
        continue;

      coords = GEOSCoordSeq_create_r(handle, shape->line[j].numpoints,
                                     2); /* todo handle z's */
      if (!coords) {
        free(innerRings);
        free(innerList);
        return NULL; /* todo, this will leak memory (shell + allocated holes) */
      }

      for (i = 0; i < shape->line[j].numpoints; i++) {
        GEOSCoordSeq_setX_r(handle, coords, i, shape->line[j].point[i].x);
        GEOSCoordSeq_setY_r(handle, coords, i, shape->line[j].point[i].y);
        /* GEOSCoordSeq_setZ(coords, i, shape->line[j].point[i].z); */
      }

      innerRings[k] = GEOSGeom_createLinearRing_r(
          handle, coords); /* innerRings[k] owns the coordinates in coords */
      k++;
    }
  }

  g = GEOSGeom_createPolygon_r(handle, outerRing, innerRings, numInnerRings);

  free(innerList);  /* clean up */
  free(innerRings); /* clean up */

  return g;
}

static GEOSGeom msGEOSShape2Geometry_polygon(shapeObj *shape) {
  int i, j;
  GEOSGeom *polygons;
  int *outerList, numOuterRings = 0, lastOuterRing = 0;
  GEOSGeom g;
  GEOSContextHandle_t handle = msGetGeosContextHandle();

  outerList = msGetOuterList(shape);
  for (i = 0; i < shape->numlines; i++) {
    if (outerList[i] == MS_TRUE) {
      numOuterRings++;
      lastOuterRing = i; /* save for the simple case */
    }
  }

  if (numOuterRings == 1) {
    g = msGEOSShape2Geometry_simplepolygon(shape, lastOuterRing, outerList);
  } else { /* a true multipolygon */
    polygons = msSmallMalloc(numOuterRings * sizeof(GEOSGeom));

    j = 0; /* part counter */
    for (i = 0; i < shape->numlines; i++) {
      if (outerList[i] == MS_FALSE)
        continue;
      polygons[j] = msGEOSShape2Geometry_simplepolygon(
          shape, i, outerList); /* TODO: account for NULL return values */
      j++;
    }

    g = GEOSGeom_createCollection_r(handle, GEOS_MULTIPOLYGON, polygons,
                                    numOuterRings);
    free(polygons);
  }

  free(outerList);
  return g;
}

GEOSGeom msGEOSShape2Geometry(shapeObj *shape) {
  if (!shape)
    return NULL; /* a NULL shape generates a NULL geometry */

  switch (shape->type) {
  case MS_SHAPE_POINT:
    if (shape->numlines == 0 ||
        shape->line[0].numpoints == 0) /* not enough info for a point */
      return NULL;

    if (shape->line[0].numpoints == 1) /* simple point */
      return msGEOSShape2Geometry_point(&(shape->line[0].point[0]));
    else /* multi-point */
      return msGEOSShape2Geometry_multipoint(&(shape->line[0]));
    break;
  case MS_SHAPE_LINE:
    if (shape->numlines == 0 ||
        shape->line[0].numpoints < 2) /* not enough info for a line */
      return NULL;

    if (shape->numlines == 1) /* simple line */
      return msGEOSShape2Geometry_line(&(shape->line[0]));
    else /* multi-line */
      return msGEOSShape2Geometry_multiline(shape);
    break;
  case MS_SHAPE_POLYGON:
    if (shape->numlines == 0 ||
        shape->line[0].numpoints <
            4) /* not enough info for a polygon (first=last) */
      return NULL;

    return msGEOSShape2Geometry_polygon(
        shape); /* simple and multipolygon cases are addressed */
    break;
  default:
    break;
  }

  return NULL; /* should not get here */
}

static shapeObj *msGEOSGeometry2Shape_point(GEOSGeom g) {
  GEOSCoordSeq coords;
  shapeObj *shape = NULL;
  GEOSContextHandle_t handle = msGetGeosContextHandle();

  if (!g)
    return NULL;

  shape = (shapeObj *)malloc(sizeof(shapeObj));
  msInitShape(shape);

  shape->type = MS_SHAPE_POINT;
  shape->line = (lineObj *)malloc(sizeof(lineObj));
  shape->numlines = 1;
  shape->line[0].point = (pointObj *)malloc(sizeof(pointObj));
  shape->line[0].numpoints = 1;
  shape->geometry = (GEOSGeom)g;

  coords = (GEOSCoordSeq)GEOSGeom_getCoordSeq_r(handle, g);

  GEOSCoordSeq_getX_r(handle, coords, 0, &(shape->line[0].point[0].x));
  GEOSCoordSeq_getY_r(handle, coords, 0, &(shape->line[0].point[0].y));
  /* GEOSCoordSeq_getZ(coords, 0, &(shape->line[0].point[0].z)); */

  shape->bounds.minx = shape->bounds.maxx = shape->line[0].point[0].x;
  shape->bounds.miny = shape->bounds.maxy = shape->line[0].point[0].y;

  return shape;
}

static shapeObj *msGEOSGeometry2Shape_multipoint(GEOSGeom g) {
  int i;
  int numPoints;
  GEOSCoordSeq coords;
  GEOSGeom point;

  shapeObj *shape = NULL;
  GEOSContextHandle_t handle = msGetGeosContextHandle();

  if (!g)
    return NULL;
  numPoints = GEOSGetNumGeometries_r(handle, g); /* each geometry has 1 point */

  shape = (shapeObj *)malloc(sizeof(shapeObj));
  msInitShape(shape);

  shape->type = MS_SHAPE_POINT;
  shape->line = (lineObj *)malloc(sizeof(lineObj));
  shape->numlines = 1;
  shape->line[0].point = (pointObj *)malloc(sizeof(pointObj) * numPoints);
  shape->line[0].numpoints = numPoints;
  shape->geometry = (GEOSGeom)g;

  for (i = 0; i < numPoints; i++) {
    point = (GEOSGeom)GEOSGetGeometryN_r(handle, g, i);
    coords = (GEOSCoordSeq)GEOSGeom_getCoordSeq_r(handle, point);

    GEOSCoordSeq_getX_r(handle, coords, 0, &(shape->line[0].point[i].x));
    GEOSCoordSeq_getY_r(handle, coords, 0, &(shape->line[0].point[i].y));
    /* GEOSCoordSeq_getZ(coords, 0, &(shape->line[0].point[i].z)); */
  }

  msComputeBounds(shape);

  return shape;
}

static shapeObj *msGEOSGeometry2Shape_line(GEOSGeom g) {
  shapeObj *shape = NULL;
  GEOSContextHandle_t handle = msGetGeosContextHandle();

  int i;
  int numPoints;
  GEOSCoordSeq coords;

  if (!g)
    return NULL;
  numPoints = GEOSGetNumCoordinates_r(handle, g);
  coords = (GEOSCoordSeq)GEOSGeom_getCoordSeq_r(handle, g);

  shape = (shapeObj *)malloc(sizeof(shapeObj));
  msInitShape(shape);

  shape->type = MS_SHAPE_LINE;
  shape->line = (lineObj *)malloc(sizeof(lineObj));
  shape->numlines = 1;
  shape->line[0].point = (pointObj *)malloc(sizeof(pointObj) * numPoints);
  shape->line[0].numpoints = numPoints;
  shape->geometry = (GEOSGeom)g;

  for (i = 0; i < numPoints; i++) {
    GEOSCoordSeq_getX_r(handle, coords, i, &(shape->line[0].point[i].x));
    GEOSCoordSeq_getY_r(handle, coords, i, &(shape->line[0].point[i].y));
    /* GEOSCoordSeq_getZ(coords, i, &(shape->line[0].point[i].z)); */
  }

  msComputeBounds(shape);

  return shape;
}

static void msGEOSGeometry2Shape_multiline_part(GEOSContextHandle_t handle,
                                                GEOSGeom part,
                                                shapeObj *shape) {
  int i;
  int type, numGeometries, numPoints;
  GEOSCoordSeq coords;
  lineObj line;

  type = GEOSGeomTypeId_r(handle, part);
  if (type == GEOS_MULTILINESTRING) {
    numGeometries = GEOSGetNumGeometries_r(handle, part);

    for (i = 0; i < numGeometries; i++) {
      GEOSGeom subPart = (GEOSGeom)GEOSGetGeometryN_r(handle, part, i);
      msGEOSGeometry2Shape_multiline_part(handle, subPart, shape);
    }
  } else {
    numPoints = GEOSGetNumCoordinates_r(handle, part);
    coords = (GEOSCoordSeq)GEOSGeom_getCoordSeq_r(handle, part);
    line.point = (pointObj *)malloc(sizeof(pointObj) * numPoints);
    line.numpoints = numPoints;

    for (i = 0; i < numPoints; i++) {
      GEOSCoordSeq_getX_r(handle, coords, i, &(line.point[i].x));
      GEOSCoordSeq_getY_r(handle, coords, i, &(line.point[i].y));
      /* GEOSCoordSeq_getZ(coords, i, &(line.point[i].z)); */
    }

    msAddLineDirectly(shape, &line);
  }
}

static shapeObj *msGEOSGeometry2Shape_multiline(GEOSGeom g) {
  int i;
  int numGeometries;
  GEOSGeom part;

  shapeObj *shape = NULL;
  GEOSContextHandle_t handle = msGetGeosContextHandle();

  if (!g)
    return NULL;
  numGeometries = GEOSGetNumGeometries_r(handle, g);

  shape = (shapeObj *)malloc(sizeof(shapeObj));
  msInitShape(shape);

  shape->type = MS_SHAPE_LINE;
  shape->geometry = (GEOSGeom)g;

  for (i = 0; i < numGeometries; i++) {
    part = (GEOSGeom)GEOSGetGeometryN_r(handle, g, i);
    msGEOSGeometry2Shape_multiline_part(handle, part, shape);
  }

  msComputeBounds(shape);

  return shape;
}

static shapeObj *msGEOSGeometry2Shape_polygon(GEOSGeom g) {
  shapeObj *shape = NULL;
  lineObj line;
  int numPoints, numRings;
  int i, j;

  GEOSCoordSeq coords;
  GEOSGeom ring;
  GEOSContextHandle_t handle = msGetGeosContextHandle();

  if (!g)
    return NULL;

  shape = (shapeObj *)malloc(sizeof(shapeObj));
  msInitShape(shape);
  shape->type = MS_SHAPE_POLYGON;
  shape->geometry = (GEOSGeom)g;

  /* exterior ring */
  ring = (GEOSGeom)GEOSGetExteriorRing_r(handle, g);
  numPoints = GEOSGetNumCoordinates_r(handle, ring);
  coords = (GEOSCoordSeq)GEOSGeom_getCoordSeq_r(handle, ring);

  line.point = (pointObj *)malloc(sizeof(pointObj) * numPoints);
  line.numpoints = numPoints;

  for (i = 0; i < numPoints; i++) {
    GEOSCoordSeq_getX_r(handle, coords, i, &(line.point[i].x));
    GEOSCoordSeq_getY_r(handle, coords, i, &(line.point[i].y));
    /* GEOSCoordSeq_getZ(coords, i, &(line.point[i].z)); */
  }
  msAddLineDirectly(shape, &line);

  /* interior rings */
  numRings = GEOSGetNumInteriorRings_r(handle, g);
  for (j = 0; j < numRings; j++) {
    ring = (GEOSGeom)GEOSGetInteriorRingN_r(handle, g, j);
    if (GEOSisRing_r(handle, ring) != 1)
      continue; /* skip it */

    numPoints = GEOSGetNumCoordinates_r(handle, ring);
    coords = (GEOSCoordSeq)GEOSGeom_getCoordSeq_r(handle, ring);

    line.point = (pointObj *)malloc(sizeof(pointObj) * numPoints);
    line.numpoints = numPoints;

    for (i = 0; i < numPoints; i++) {
      GEOSCoordSeq_getX_r(handle, coords, i, &(line.point[i].x));
      GEOSCoordSeq_getY_r(handle, coords, i, &(line.point[i].y));
      /* GEOSCoordSeq_getZ(coords, i, &(line.point[i].z)); */
    }
    msAddLineDirectly(shape, &line);
  }

  msComputeBounds(shape);

  return shape;
}

static shapeObj *msGEOSGeometry2Shape_multipolygon(GEOSGeom g) {
  int i, j, k;
  shapeObj *shape = NULL;
  lineObj line;
  int numPoints, numRings, numPolygons;

  GEOSCoordSeq coords;
  GEOSGeom polygon, ring;
  GEOSContextHandle_t handle = msGetGeosContextHandle();

  if (!g)
    return NULL;
  numPolygons = GEOSGetNumGeometries_r(handle, g);

  shape = (shapeObj *)malloc(sizeof(shapeObj));
  msInitShape(shape);
  shape->type = MS_SHAPE_POLYGON;
  shape->geometry = (GEOSGeom)g;

  for (k = 0; k < numPolygons; k++) { /* for each polygon */
    polygon = (GEOSGeom)GEOSGetGeometryN_r(handle, g, k);

    /* exterior ring */
    ring = (GEOSGeom)GEOSGetExteriorRing_r(handle, polygon);
    numPoints = GEOSGetNumCoordinates_r(handle, ring);
    coords = (GEOSCoordSeq)GEOSGeom_getCoordSeq_r(handle, ring);

    line.point = (pointObj *)malloc(sizeof(pointObj) * numPoints);
    line.numpoints = numPoints;

    for (i = 0; i < numPoints; i++) {
      GEOSCoordSeq_getX_r(handle, coords, i, &(line.point[i].x));
      GEOSCoordSeq_getY_r(handle, coords, i, &(line.point[i].y));
      /* GEOSCoordSeq_getZ(coords, i, &(line.point[i].z)); */
    }
    msAddLineDirectly(shape, &line);

    /* interior rings */
    numRings = GEOSGetNumInteriorRings_r(handle, polygon);

    for (j = 0; j < numRings; j++) {
      ring = (GEOSGeom)GEOSGetInteriorRingN_r(handle, polygon, j);
      if (GEOSisRing_r(handle, ring) != 1)
        continue; /* skip it */

      numPoints = GEOSGetNumCoordinates_r(handle, ring);
      coords = (GEOSCoordSeq)GEOSGeom_getCoordSeq_r(handle, ring);

      line.point = (pointObj *)malloc(sizeof(pointObj) * numPoints);
      line.numpoints = numPoints;

      for (i = 0; i < numPoints; i++) {
        GEOSCoordSeq_getX_r(handle, coords, i, &(line.point[i].x));
        GEOSCoordSeq_getY_r(handle, coords, i, &(line.point[i].y));
        /* GEOSCoordSeq_getZ(coords, i, &(line.point[i].z)); */
      }
      msAddLineDirectly(shape, &line);
    }
  } /* next polygon */

  msComputeBounds(shape);

  return shape;
}

shapeObj *msGEOSGeometry2Shape(GEOSGeom g) {
  int type;
  GEOSContextHandle_t handle = msGetGeosContextHandle();

  if (!g)
    return NULL; /* a NULL geometry generates a NULL shape */

  type = GEOSGeomTypeId_r(handle, g);
  switch (type) {
  case GEOS_POINT:
    return msGEOSGeometry2Shape_point(g);
    break;
  case GEOS_MULTIPOINT:
    return msGEOSGeometry2Shape_multipoint(g);
    break;
  case GEOS_LINESTRING:
    return msGEOSGeometry2Shape_line(g);
    break;
  case GEOS_MULTILINESTRING:
    return msGEOSGeometry2Shape_multiline(g);
    break;
  case GEOS_POLYGON:
    return msGEOSGeometry2Shape_polygon(g);
    break;
  case GEOS_MULTIPOLYGON:
    return msGEOSGeometry2Shape_multipolygon(g);
    break;
  case GEOS_GEOMETRYCOLLECTION:
    if (!GEOSisEmpty_r(handle, g)) {
      shapeObj *shape = (shapeObj *)malloc(sizeof(shapeObj));
      msInitShape(shape);
      shape->type = MS_SHAPE_LINE;
      shape->geometry = (GEOSGeom)g;

      const int numGeoms = GEOSGetNumGeometries_r(handle, g);
      for (int i = 0; i < numGeoms; i++) { /* for each geometry */
        shapeObj *shape2 =
            msGEOSGeometry2Shape((GEOSGeom)GEOSGetGeometryN_r(handle, g, i));
        if (shape2) {
          for (int j = 0; j < shape2->numlines; j++)
            msAddLineDirectly(shape, &shape2->line[j]);
          shape2->numlines = 0;
          shape2->geometry = NULL; /* not owned */
          msFreeShape(shape2);
        }
      }
      msComputeBounds(shape);
      return shape;
    }
    break;
  default:
    msSetError(MS_GEOSERR, "Unsupported GEOS geometry type (%d).",
               "msGEOSGeometry2Shape()", type);
  }
  return NULL;
}
#endif

/*
** Maintenence functions exposed to MapServer/MapScript.
*/

void msGEOSFreeGeometry(shapeObj *shape) {
#ifdef USE_GEOS
  GEOSGeom g = NULL;
  GEOSContextHandle_t handle = msGetGeosContextHandle();

  if (!shape || !shape->geometry)
    return;

  g = (GEOSGeom)shape->geometry;
  GEOSGeom_destroy_r(handle, g);
  shape->geometry = NULL;
#else
  msSetError(MS_GEOSERR, "GEOS support is not available.",
             "msGEOSFreeGEOSGeom()");
  return;
#endif
}

/*
** WKT input and output functions
*/
shapeObj *msGEOSShapeFromWKT(const char *wkt) {
#ifdef USE_GEOS
  GEOSGeom g;
  GEOSContextHandle_t handle = msGetGeosContextHandle();

  if (!wkt)
    return NULL;

  g = GEOSGeomFromWKT_r(handle, wkt);
  if (!g) {
    msSetError(MS_GEOSERR, "Error reading WKT geometry \"%s\".",
               "msGEOSShapeFromWKT()", wkt);
    return NULL;
  } else {
    return msGEOSGeometry2Shape(g);
  }

#else
  msSetError(MS_GEOSERR, "GEOS support is not available.",
             "msGEOSShapeFromWKT()");
  return NULL;
#endif
}

/* Return should be freed with msGEOSFreeWKT */
char *msGEOSShapeToWKT(shapeObj *shape) {
#ifdef USE_GEOS
  GEOSGeom g;
  GEOSContextHandle_t handle = msGetGeosContextHandle();

  if (!shape)
    return NULL;

  /* if we have a geometry, we should update it*/
  msGEOSFreeGeometry(shape);

  shape->geometry = (GEOSGeom)msGEOSShape2Geometry(shape);
  g = (GEOSGeom)shape->geometry;
  if (!g)
    return NULL;

  return GEOSGeomToWKT_r(handle, g);
#else
  msSetError(MS_GEOSERR, "GEOS support is not available.",
             "msGEOSShapeToWKT()");
  return NULL;
#endif
}

void msGEOSFreeWKT(char *pszGEOSWKT) {
#ifdef USE_GEOS
  GEOSContextHandle_t handle = msGetGeosContextHandle();
#if GEOS_VERSION_MAJOR > 3 ||                                                  \
    (GEOS_VERSION_MAJOR == 3 && GEOS_VERSION_MINOR >= 2)
  GEOSFree_r(handle, pszGEOSWKT);
#endif
#else
  msSetError(MS_GEOSERR, "GEOS support is not available.", "msGEOSFreeWKT()");
#endif
}

shapeObj *msGEOSOffsetCurve(shapeObj *p, double offset) {
#if defined USE_GEOS && (GEOS_VERSION_MAJOR > 3 ||                             \
                         (GEOS_VERSION_MAJOR == 3 && GEOS_VERSION_MINOR >= 3))
  int typeChanged = 0;
  GEOSGeom g1, g2 = NULL;
  GEOSContextHandle_t handle = msGetGeosContextHandle();

  if (!p)
    return NULL;

  /*
   * GEOSOffsetCurve_r() uses BufferBuilder.bufferLineSingleSided(), which
   * works with lines, naturally. In order to allow offsets for a MapServer
   * polygonObj, it has to be processed as line and afterwards reverted.
   */
  if (p->type == MS_SHAPE_POLYGON) {
    p->type = MS_SHAPE_LINE;
    typeChanged = 1;
    msGEOSFreeGeometry(p);
  }

  if (typeChanged || !p->geometry)
    p->geometry = (GEOSGeom)msGEOSShape2Geometry(p);

  g1 = (GEOSGeom)p->geometry;
  if (!g1)
    return NULL;

  if (GEOSGeomTypeId_r(handle, g1) == GEOS_MULTILINESTRING) {
    GEOSGeom *lines = malloc(p->numlines * sizeof(GEOSGeom));
    if (!lines)
      return NULL;
    for (int i = 0; i < p->numlines; i++) {
      lines[i] =
          GEOSOffsetCurve_r(handle, GEOSGetGeometryN_r(handle, g1, i), offset,
                            4, GEOSBUF_JOIN_MITRE, fabs(offset * 1.5));
    }
    g2 = GEOSGeom_createCollection_r(handle, GEOS_MULTILINESTRING, lines,
                                     p->numlines);
    free(lines);
  } else {
    g2 = GEOSOffsetCurve_r(handle, g1, offset, 4, GEOSBUF_JOIN_MITRE,
                           fabs(offset * 1.5));
  }

  /*
   * Undo change of geometry type. We won't re-create the geos gemotry here,
   * it's up to each geos function to create it.
   */
  if (typeChanged) {
    msGEOSFreeGeometry(p);
    p->type = MS_SHAPE_POLYGON;
  }

  if (g2)
    return msGEOSGeometry2Shape(g2);

  return NULL;
#else
  msSetError(MS_GEOSERR, "GEOS Offset Curve support is not available.",
             "msGEOSingleSidedBuffer()");
  return NULL;
#endif
}

/*
** Analytical functions exposed to MapServer/MapScript.
*/

shapeObj *msGEOSBuffer(shapeObj *shape, double width) {
#ifdef USE_GEOS
  GEOSGeom g1, g2;
  GEOSContextHandle_t handle = msGetGeosContextHandle();

  if (!shape)
    return NULL;

  if (!shape->geometry) /* if no geometry for the shape then build one */
    shape->geometry = (GEOSGeom)msGEOSShape2Geometry(shape);

  g1 = (GEOSGeom)shape->geometry;
  if (!g1)
    return NULL;

  g2 = GEOSBuffer_r(handle, g1, width, 30);
  return msGEOSGeometry2Shape(g2);
#else
  msSetError(MS_GEOSERR, "GEOS support is not available.", "msGEOSBuffer()");
  return NULL;
#endif
}

shapeObj *msGEOSSimplify(shapeObj *shape, double tolerance) {
#ifdef USE_GEOS
  GEOSGeom g1, g2;
  GEOSContextHandle_t handle = msGetGeosContextHandle();

  if (!shape)
    return NULL;

  if (!shape->geometry) /* if no geometry for the shape then build one */
    shape->geometry = (GEOSGeom)msGEOSShape2Geometry(shape);

  g1 = (GEOSGeom)shape->geometry;
  if (!g1)
    return NULL;

  g2 = GEOSSimplify_r(handle, g1, tolerance);
  return msGEOSGeometry2Shape(g2);
#else
  msSetError(MS_GEOSERR, "GEOS Simplify support is not available.",
             "msGEOSSimplify()");
  return NULL;
#endif
}

shapeObj *msGEOSTopologyPreservingSimplify(shapeObj *shape, double tolerance) {
#ifdef USE_GEOS
  GEOSGeom g1, g2;
  GEOSContextHandle_t handle = msGetGeosContextHandle();

  if (!shape)
    return NULL;

  if (!shape->geometry) /* if no geometry for the shape then build one */
    shape->geometry = (GEOSGeom)msGEOSShape2Geometry(shape);

  g1 = (GEOSGeom)shape->geometry;
  if (!g1)
    return NULL;

  g2 = GEOSTopologyPreserveSimplify_r(handle, g1, tolerance);
  return msGEOSGeometry2Shape(g2);
#else
  msSetError(MS_GEOSERR, "GEOS Simplify support is not available.",
             "msGEOSTopologyPreservingSimplify()");
  return NULL;
#endif
}

shapeObj *msGEOSConvexHull(shapeObj *shape) {
#ifdef USE_GEOS
  GEOSGeom g1, g2;
  GEOSContextHandle_t handle = msGetGeosContextHandle();

  if (!shape)
    return NULL;

  if (!shape->geometry) /* if no geometry for the shape then build one */
    shape->geometry = (GEOSGeom)msGEOSShape2Geometry(shape);
  g1 = (GEOSGeom)shape->geometry;
  if (!g1)
    return NULL;

  g2 = GEOSConvexHull_r(handle, g1);
  return msGEOSGeometry2Shape(g2);
#else
  msSetError(MS_GEOSERR, "GEOS support is not available.",
             "msGEOSConvexHull()");
  return NULL;
#endif
}

shapeObj *msGEOSBoundary(shapeObj *shape) {
#ifdef USE_GEOS
  GEOSGeom g1, g2;
  GEOSContextHandle_t handle = msGetGeosContextHandle();

  if (!shape)
    return NULL;

  if (!shape->geometry) /* if no geometry for the shape then build one */
    shape->geometry = (GEOSGeom)msGEOSShape2Geometry(shape);
  g1 = (GEOSGeom)shape->geometry;
  if (!g1)
    return NULL;

  g2 = GEOSBoundary_r(handle, g1);
  return msGEOSGeometry2Shape(g2);
#else
  msSetError(MS_GEOSERR, "GEOS support is not available.", "msGEOSBoundary()");
  return NULL;
#endif
}

pointObj *msGEOSGetCentroid(shapeObj *shape) {
#ifdef USE_GEOS
  GEOSGeom g1, g2;
  GEOSCoordSeq coords;
  pointObj *point;
  GEOSContextHandle_t handle = msGetGeosContextHandle();

  if (!shape)
    return NULL;

  if (!shape->geometry) /* if no geometry for the shape then build one */
    shape->geometry = (GEOSGeom)msGEOSShape2Geometry(shape);
  g1 = (GEOSGeom)shape->geometry;
  if (!g1)
    return NULL;

  g2 = GEOSGetCentroid_r(handle, g1);
  if (!g2)
    return NULL;

  point = (pointObj *)malloc(sizeof(pointObj));

  coords = (GEOSCoordSeq)GEOSGeom_getCoordSeq_r(handle, g2);

  GEOSCoordSeq_getX_r(handle, coords, 0, &(point->x));
  GEOSCoordSeq_getY_r(handle, coords, 0, &(point->y));
  /* GEOSCoordSeq_getZ(coords, 0, &(point->z)); */

  GEOSGeom_destroy_r(handle, g2);

  return point;
#else
  msSetError(MS_GEOSERR, "GEOS support is not available.",
             "msGEOSGetCentroid()");
  return NULL;
#endif
}

shapeObj *msGEOSUnion(shapeObj *shape1, shapeObj *shape2) {
#ifdef USE_GEOS
  GEOSGeom g1, g2, g3;
  GEOSContextHandle_t handle = msGetGeosContextHandle();

  if (!shape1 || !shape2)
    return NULL;

  if (!shape1->geometry) /* if no geometry for the shape then build one */
    shape1->geometry = (GEOSGeom)msGEOSShape2Geometry(shape1);
  g1 = (GEOSGeom)shape1->geometry;
  if (!g1)
    return NULL;

  if (!shape2->geometry) /* if no geometry for the shape then build one */
    shape2->geometry = (GEOSGeom)msGEOSShape2Geometry(shape2);
  g2 = (GEOSGeom)shape2->geometry;
  if (!g2)
    return NULL;

  g3 = GEOSUnion_r(handle, g1, g2);
  return msGEOSGeometry2Shape(g3);
#else
  msSetError(MS_GEOSERR, "GEOS support is not available.", "msGEOSUnion()");
  return NULL;
#endif
}

shapeObj *msGEOSIntersection(shapeObj *shape1, shapeObj *shape2) {
#ifdef USE_GEOS
  GEOSGeom g1, g2, g3;
  GEOSContextHandle_t handle = msGetGeosContextHandle();

  if (!shape1 || !shape2)
    return NULL;

  if (!shape1->geometry) /* if no geometry for the shape then build one */
    shape1->geometry = (GEOSGeom)msGEOSShape2Geometry(shape1);
  g1 = (GEOSGeom)shape1->geometry;
  if (!g1)
    return NULL;

  if (!shape2->geometry) /* if no geometry for the shape then build one */
    shape2->geometry = (GEOSGeom)msGEOSShape2Geometry(shape2);
  g2 = (GEOSGeom)shape2->geometry;
  if (!g2)
    return NULL;

  g3 = GEOSIntersection_r(handle, g1, g2);
  return msGEOSGeometry2Shape(g3);
#else
  msSetError(MS_GEOSERR, "GEOS support is not available.",
             "msGEOSIntersection()");
  return NULL;
#endif
}

shapeObj *msGEOSDifference(shapeObj *shape1, shapeObj *shape2) {
#ifdef USE_GEOS
  GEOSGeom g1, g2, g3;
  GEOSContextHandle_t handle = msGetGeosContextHandle();

  if (!shape1 || !shape2)
    return NULL;

  if (!shape1->geometry) /* if no geometry for the shape then build one */
    shape1->geometry = (GEOSGeom)msGEOSShape2Geometry(shape1);
  g1 = (GEOSGeom)shape1->geometry;
  if (!g1)
    return NULL;

  if (!shape2->geometry) /* if no geometry for the shape then build one */
    shape2->geometry = (GEOSGeom)msGEOSShape2Geometry(shape2);
  g2 = (GEOSGeom)shape2->geometry;
  if (!g2)
    return NULL;

  g3 = GEOSDifference_r(handle, g1, g2);
  return msGEOSGeometry2Shape(g3);
#else
  msSetError(MS_GEOSERR, "GEOS support is not available.",
             "msGEOSDifference()");
  return NULL;
#endif
}

shapeObj *msGEOSSymDifference(shapeObj *shape1, shapeObj *shape2) {
#ifdef USE_GEOS
  GEOSGeom g1, g2, g3;
  GEOSContextHandle_t handle = msGetGeosContextHandle();

  if (!shape1 || !shape2)
    return NULL;

  if (!shape1->geometry) /* if no geometry for the shape then build one */
    shape1->geometry = (GEOSGeom)msGEOSShape2Geometry(shape1);
  g1 = (GEOSGeom)shape1->geometry;
  if (!g1)
    return NULL;

  if (!shape2->geometry) /* if no geometry for the shape then build one */
    shape2->geometry = (GEOSGeom)msGEOSShape2Geometry(shape2);
  g2 = (GEOSGeom)shape2->geometry;
  if (!g2)
    return NULL;

  g3 = GEOSSymDifference_r(handle, g1, g2);
  return msGEOSGeometry2Shape(g3);
#else
  msSetError(MS_GEOSERR, "GEOS support is not available.",
             "msGEOSSymDifference()");
  return NULL;
#endif
}

shapeObj *msGEOSLineMerge(shapeObj *shape) {
#ifdef USE_GEOS
  GEOSGeom g1, g2;
  GEOSContextHandle_t handle = msGetGeosContextHandle();

  if (!shape)
    return NULL;
  if (shape->type != MS_SHAPE_LINE)
    return NULL;

  if (!shape->geometry) /* if no geometry for the shape then build one */
    shape->geometry = (GEOSGeom)msGEOSShape2Geometry(shape);
  g1 = (GEOSGeom)shape->geometry;
  if (!g1)
    return NULL;

  g2 = GEOSLineMerge_r(handle, g1);
  return msGEOSGeometry2Shape(g2);
#else
  msSetError(MS_GEOSERR, "GEOS support is not available.", "msGEOSLineMerge()");
  return NULL;
#endif
}

shapeObj *msGEOSVoronoiDiagram(shapeObj *shape, double tolerance,
                               int onlyEdges) {
#if defined(USE_GEOS) &&                                                       \
    (GEOS_VERSION_MAJOR > 3 ||                                                 \
     (GEOS_VERSION_MAJOR == 3 && GEOS_VERSION_MINOR >= 5))
  GEOSGeom g1, g2;
  GEOSContextHandle_t handle = msGetGeosContextHandle();

  if (!shape)
    return NULL;

  if (!shape->geometry) /* if no geometry for the shape then build one */
    shape->geometry = (GEOSGeom)msGEOSShape2Geometry(shape);
  g1 = (GEOSGeom)shape->geometry;
  if (!g1)
    return NULL;

  g2 = GEOSVoronoiDiagram_r(handle, g1, NULL, tolerance, onlyEdges);
  return msGEOSGeometry2Shape(g2);
#else
  msSetError(
      MS_GEOSERR,
      "GEOS support is not available or GEOS version is not 3.5 or higher.",
      "msGEOSVoronoiDiagram()");
  return NULL;
#endif
}

static int keepEdge(lineObj *segment, shapeObj *polygon) {
  int i, j;

  if (segment->numpoints < 2)
    return MS_FALSE;
  if (msIntersectPointPolygon(&segment->point[0], polygon) != MS_TRUE)
    return MS_FALSE;
  if (msIntersectPointPolygon(&segment->point[1], polygon) != MS_TRUE)
    return MS_FALSE;

  for (i = 0; i < polygon->numlines; i++)
    for (j = 1; j < polygon->line[i].numpoints; j++)
      if (msIntersectSegments(&(segment->point[0]), &(segment->point[1]),
                              &(polygon->line[i].point[j - 1]),
                              &(polygon->line[i].point[j])) == MS_TRUE)
        return (MS_FALSE);

  return (MS_TRUE);
}

#define ARE_SAME_POINTS(a, b)                                                  \
  (((a).x != (b).x || (a).y != (b).y) ? MS_FALSE : MS_TRUE)

// returns the index of the node, we use z to store a count of points at the
// same coordinate
static int buildNodes(multipointObj *nodes, pointObj *point) {
  int i;

  for (i = 0; i < nodes->numpoints; i++) {
    if (ARE_SAME_POINTS(nodes->point[i], *point)) { // found it
      nodes->point[i].z++;
      return i;
    }
  }

  // not found, add it
  nodes->point[i].x = point->x;
  nodes->point[i].y = point->y;
  nodes->point[i].z = 1;
  nodes->numpoints++;

  return i;
}

shapeObj *msGEOSCenterline(shapeObj *shape) {
#if defined(USE_GEOS) &&                                                       \
    (GEOS_VERSION_MAJOR > 3 ||                                                 \
     (GEOS_VERSION_MAJOR == 3 && GEOS_VERSION_MINOR >= 5))
  int i;
  shapeObj *shape2 = NULL;

  multipointObj nodes;
  graphObj *graph;

  int *path = NULL; // array of node indexes
  int path_size = 0;
  double path_dist = -1, max_path_dist = -1;

  if (!shape)
    return NULL;
  if (shape->type != MS_SHAPE_POLYGON) {
    msSetError(MS_GEOSERR,
               "Centerlines can only be computed for polygon shapes.",
               "msGEOSCenterline()");
    return NULL;
  }

  shape2 = msGEOSVoronoiDiagram(shape, 0.0, MS_TRUE);
  if (!shape2) {
    msSetError(MS_GEOSERR, "Voronoi diagram generation failed.",
               "msGEOSCenterline()");
    return NULL;
  }

  // process the edges and build a graph representation
  nodes.point = (pointObj *)malloc(shape2->numlines * sizeof(pointObj) * 2);
  nodes.numpoints = 0;
  if (!nodes.point) {
    msFreeShape(shape2);
    free(shape2);
    return NULL;
  }

  graph = msCreateGraph(shape2->numlines * 2);
  if (!graph) {
    msFreeShape(shape2);
    free(shape2);
    free(nodes.point);
    return NULL;
  }

  for (i = 0; i < shape2->numlines; i++) {
    if (keepEdge(&shape2->line[i], shape) == MS_TRUE) {
      int src = buildNodes(&nodes, &shape2->line[i].point[0]);
      int dest = buildNodes(&nodes, &shape2->line[i].point[1]);
      msGraphAddEdge(graph, src, dest,
                     msDistancePointToPoint(&shape2->line[i].point[0],
                                            &shape2->line[i].point[1]));
    }
  }
  msFreeShape(shape2); // done with voronoi geometry, shape2 is still allocated
                       // though, just empty
  shape2->type = MS_SHAPE_LINE; // will fill with centerline

  if (nodes.numpoints == 0) {
    msSetError(MS_GEOSERR,
               "Centerline generation failed, try densifying the shapes.",
               "msGEOSCenterline()");
    free(shape2);
    msFreeGraph(graph);
    free(nodes.point);
    return NULL;
  }

  // step through edge nodes (z=1)
  for (i = 0; i < nodes.numpoints; i++) {
    if (nodes.point[i].z != 1)
      continue; // skip

    if (!path) { // first one, keep this path
      path = msGraphGetLongestShortestPath(graph, i, &path_size, &path_dist);
      max_path_dist = path_dist;
    } else {
      int *tmp_path = NULL;
      int tmp_path_size = 0;
      double tmp_path_dist = -1;

      if (i == path[path_size - 1])
        continue; // skip, graph is bi-directional so it can't be any longer
      tmp_path = msGraphGetLongestShortestPath(graph, i, &tmp_path_size,
                                               &tmp_path_dist);
      if (tmp_path_dist > max_path_dist) {
        free(path);
        path = tmp_path;
        path_size = tmp_path_size;
        path_dist = tmp_path_dist;
        max_path_dist = tmp_path_dist;
      } else { // skip path
        msFree(tmp_path);
      }
    }
  }
  msFreeGraph(graph); // done with graph

  // transform the path into a shape
  if (!path) {
    msSetError(MS_GEOSERR, "Centerline generation failed.",
               "msGEOSCenterline()");
    free(shape2);
    free(nodes.point);
    return NULL;
  } else {
    lineObj line;
    line.point = (pointObj *)malloc(path_size * sizeof(pointObj));
    if (!line.point) {
      free(shape2);
      free(path);
      free(nodes.point);
      return NULL;
    }
    line.numpoints = path_size;

    for (i = 0; i < path_size; i++) {
      line.point[i].x = nodes.point[path[i]].x;
      line.point[i].y = nodes.point[path[i]].y;
    }
    msAddLineDirectly(shape2, &line);
  }

  free(path); // clean up
  free(nodes.point);

  return shape2;
#else
  msSetError(
      MS_GEOSERR,
      "GEOS support is not available or GEOS version is not 3.5 or higher.",
      "msGEOSCenterline()");
  return NULL;
#endif
}

/*
** Binary predicates exposed to MapServer/MapScript
*/

/*
** Does shape1 contain shape2, returns MS_TRUE/MS_FALSE or -1 for an error.
*/
int msGEOSContains(shapeObj *shape1, shapeObj *shape2) {
#ifdef USE_GEOS
  GEOSGeom g1, g2;
  int result;
  GEOSContextHandle_t handle = msGetGeosContextHandle();

  if (!shape1 || !shape2)
    return -1;

  if (!shape1->geometry) /* if no geometry for shape1 then build one */
    shape1->geometry = (GEOSGeom)msGEOSShape2Geometry(shape1);
  g1 = shape1->geometry;
  if (!g1)
    return -1;

  if (!shape2->geometry) /* if no geometry for shape2 then build one */
    shape2->geometry = (GEOSGeom)msGEOSShape2Geometry(shape2);
  g2 = shape2->geometry;
  if (!g2)
    return -1;

  result = GEOSContains_r(handle, g1, g2);
  return ((result == 2) ? -1 : result);
#else
  msSetError(MS_GEOSERR, "GEOS support is not available.", "msGEOSContains()");
  return -1;
#endif
}

/*
** Does shape1 overlap shape2, returns MS_TRUE/MS_FALSE or -1 for an error.
*/
int msGEOSOverlaps(shapeObj *shape1, shapeObj *shape2) {
#ifdef USE_GEOS
  GEOSGeom g1, g2;
  int result;
  GEOSContextHandle_t handle = msGetGeosContextHandle();

  if (!shape1 || !shape2)
    return -1;

  if (!shape1->geometry) /* if no geometry for shape1 then build one */
    shape1->geometry = (GEOSGeom)msGEOSShape2Geometry(shape1);
  g1 = shape1->geometry;
  if (!g1)
    return -1;

  if (!shape2->geometry) /* if no geometry for shape2 then build one */
    shape2->geometry = (GEOSGeom)msGEOSShape2Geometry(shape2);
  g2 = shape2->geometry;
  if (!g2)
    return -1;

  result = GEOSOverlaps_r(handle, g1, g2);
  return ((result == 2) ? -1 : result);
#else
  msSetError(MS_GEOSERR, "GEOS support is not available.", "msGEOSOverlaps()");
  return -1;
#endif
}

/*
** Is shape1 within shape2, returns MS_TRUE/MS_FALSE or -1 for an error.
*/
int msGEOSWithin(shapeObj *shape1, shapeObj *shape2) {
#ifdef USE_GEOS
  GEOSGeom g1, g2;
  int result;
  GEOSContextHandle_t handle = msGetGeosContextHandle();

  if (!shape1 || !shape2)
    return -1;

  if (!shape1->geometry) /* if no geometry for shape1 then build one */
    shape1->geometry = (GEOSGeom)msGEOSShape2Geometry(shape1);
  g1 = shape1->geometry;
  if (!g1)
    return -1;

  if (!shape2->geometry) /* if no geometry for shape2 then build one */
    shape2->geometry = (GEOSGeom)msGEOSShape2Geometry(shape2);
  g2 = shape2->geometry;
  if (!g2)
    return -1;

  result = GEOSWithin_r(handle, g1, g2);
  return ((result == 2) ? -1 : result);
#else
  msSetError(MS_GEOSERR, "GEOS support is not available.", "msGEOSWithin()");
  return -1;
#endif
}

/*
** Does shape1 cross shape2, returns MS_TRUE/MS_FALSE or -1 for an error.
*/
int msGEOSCrosses(shapeObj *shape1, shapeObj *shape2) {
#ifdef USE_GEOS
  GEOSGeom g1, g2;
  int result;
  GEOSContextHandle_t handle = msGetGeosContextHandle();

  if (!shape1 || !shape2)
    return -1;

  if (!shape1->geometry) /* if no geometry for shape1 then build one */
    shape1->geometry = (GEOSGeom)msGEOSShape2Geometry(shape1);
  g1 = shape1->geometry;
  if (!g1)
    return -1;

  if (!shape2->geometry) /* if no geometry for shape2 then build one */
    shape2->geometry = (GEOSGeom)msGEOSShape2Geometry(shape2);
  g2 = shape2->geometry;
  if (!g2)
    return -1;

  result = GEOSCrosses_r(handle, g1, g2);
  return ((result == 2) ? -1 : result);
#else
  msSetError(MS_GEOSERR, "GEOS support is not available.", "msGEOSCrosses()");
  return -1;
#endif
}

/*
** Does shape1 intersect shape2, returns MS_TRUE/MS_FALSE or -1 for an error.
*/
int msGEOSIntersects(shapeObj *shape1, shapeObj *shape2) {
#ifdef USE_GEOS
  GEOSGeom g1, g2;
  int result;
  GEOSContextHandle_t handle = msGetGeosContextHandle();

  if (!shape1 || !shape2)
    return -1;

  if (!shape1->geometry) /* if no geometry for shape1 then build one */
    shape1->geometry = (GEOSGeom)msGEOSShape2Geometry(shape1);
  g1 = (GEOSGeom)shape1->geometry;
  if (!g1)
    return -1;

  if (!shape2->geometry) /* if no geometry for shape2 then build one */
    shape2->geometry = (GEOSGeom)msGEOSShape2Geometry(shape2);
  g2 = (GEOSGeom)shape2->geometry;
  if (!g2)
    return -1;

  result = GEOSIntersects_r(handle, g1, g2);
  return ((result == 2) ? -1 : result);
#else
  if (!shape1 || !shape2)
    return -1;

  switch (shape1->type) { /* todo: deal with point shapes */
  case (MS_SHAPE_LINE):
    switch (shape2->type) {
    case (MS_SHAPE_LINE):
      return msIntersectPolylines(shape1, shape2);
    case (MS_SHAPE_POLYGON):
      return msIntersectPolylinePolygon(shape1, shape2);
    }
    break;
  case (MS_SHAPE_POLYGON):
    switch (shape2->type) {
    case (MS_SHAPE_LINE):
      return msIntersectPolylinePolygon(shape2, shape1);
    case (MS_SHAPE_POLYGON):
      return msIntersectPolygons(shape1, shape2);
    }
    break;
  }

  return -1;
#endif
}

/*
** Does shape1 touch shape2, returns MS_TRUE/MS_FALSE or -1 for an error.
*/
int msGEOSTouches(shapeObj *shape1, shapeObj *shape2) {
#ifdef USE_GEOS
  GEOSGeom g1, g2;
  int result;
  GEOSContextHandle_t handle = msGetGeosContextHandle();

  if (!shape1 || !shape2)
    return -1;

  if (!shape1->geometry) /* if no geometry for shape1 then build one */
    shape1->geometry = (GEOSGeom)msGEOSShape2Geometry(shape1);
  g1 = (GEOSGeom)shape1->geometry;
  if (!g1)
    return -1;

  if (!shape2->geometry) /* if no geometry for shape2 then build one */
    shape2->geometry = (GEOSGeom)msGEOSShape2Geometry(shape2);
  g2 = (GEOSGeom)shape2->geometry;
  if (!g2)
    return -1;

  result = GEOSTouches_r(handle, g1, g2);
  return ((result == 2) ? -1 : result);
#else
  msSetError(MS_GEOSERR, "GEOS support is not available.", "msGEOSTouches()");
  return -1;
#endif
}

/*
** Does shape1 equal shape2, returns MS_TRUE/MS_FALSE or -1 for an error.
*/
int msGEOSEquals(shapeObj *shape1, shapeObj *shape2) {
#ifdef USE_GEOS
  GEOSGeom g1, g2;
  int result;
  GEOSContextHandle_t handle = msGetGeosContextHandle();

  if (!shape1 || !shape2)
    return -1;

  if (!shape1->geometry) /* if no geometry for shape1 then build one */
    shape1->geometry = (GEOSGeom)msGEOSShape2Geometry(shape1);
  g1 = (GEOSGeom)shape1->geometry;
  if (!g1)
    return -1;

  if (!shape2->geometry) /* if no geometry for shape2 then build one */
    shape2->geometry = (GEOSGeom)msGEOSShape2Geometry(shape2);
  g2 = (GEOSGeom)shape2->geometry;
  if (!g2)
    return -1;

  result = GEOSEquals_r(handle, g1, g2);
  return ((result == 2) ? -1 : result);
#else
  msSetError(MS_GEOSERR, "GEOS support is not available.", "msGEOSEquals()");
  return -1;
#endif
}

/*
** Are shape1 and shape2 disjoint, returns MS_TRUE/MS_FALSE or -1 for an error.
*/
int msGEOSDisjoint(shapeObj *shape1, shapeObj *shape2) {
#ifdef USE_GEOS
  GEOSGeom g1, g2;
  int result;
  GEOSContextHandle_t handle = msGetGeosContextHandle();

  if (!shape1 || !shape2)
    return -1;

  if (!shape1->geometry) /* if no geometry for shape1 then build one */
    shape1->geometry = (GEOSGeom)msGEOSShape2Geometry(shape1);
  g1 = (GEOSGeom)shape1->geometry;
  if (!g1)
    return -1;

  if (!shape2->geometry) /* if no geometry for shape2 then build one */
    shape2->geometry = (GEOSGeom)msGEOSShape2Geometry(shape2);
  g2 = (GEOSGeom)shape2->geometry;
  if (!g2)
    return -1;

  result = GEOSDisjoint_r(handle, g1, g2);
  return ((result == 2) ? -1 : result);
#else
  msSetError(MS_GEOSERR, "GEOS support is not available.", "msGEOSDisjoint()");
  return -1;
#endif
}

/*
** Useful misc. functions that return -1 on failure.
*/
double msGEOSArea(shapeObj *shape) {
#if defined(USE_GEOS) && defined(GEOS_CAPI_VERSION_MAJOR) &&                   \
    defined(GEOS_CAPI_VERSION_MINOR) &&                                        \
    (GEOS_CAPI_VERSION_MAJOR > 1 || GEOS_CAPI_VERSION_MINOR >= 1)
  GEOSGeom g;
  double area;
  int result;
  GEOSContextHandle_t handle = msGetGeosContextHandle();

  if (!shape)
    return -1;

  if (!shape->geometry) /* if no geometry for the shape then build one */
    shape->geometry = (GEOSGeom)msGEOSShape2Geometry(shape);
  g = (GEOSGeom)shape->geometry;
  if (!g)
    return -1;

  result = GEOSArea_r(handle, g, &area);
  return ((result == 0) ? -1 : area);
#elif defined(USE_GEOS)
  msSetError(MS_GEOSERR,
             "GEOS support enabled, but old version lacks GEOSArea().",
             "msGEOSArea()");
  return -1;
#else
  msSetError(MS_GEOSERR, "GEOS support is not available.", "msGEOSArea()");
  return -1;
#endif
}

double msGEOSLength(shapeObj *shape) {
#if defined(USE_GEOS) && defined(GEOS_CAPI_VERSION_MAJOR) &&                   \
    defined(GEOS_CAPI_VERSION_MINOR) &&                                        \
    (GEOS_CAPI_VERSION_MAJOR > 1 || GEOS_CAPI_VERSION_MINOR >= 1)

  GEOSGeom g;
  double length;
  int result;
  GEOSContextHandle_t handle = msGetGeosContextHandle();

  if (!shape)
    return -1;

  if (!shape->geometry) /* if no geometry for the shape then build one */
    shape->geometry = (GEOSGeom)msGEOSShape2Geometry(shape);
  g = (GEOSGeom)shape->geometry;
  if (!g)
    return -1;

  result = GEOSLength_r(handle, g, &length);
  return ((result == 0) ? -1 : length);
#elif defined(USE_GEOS)
  msSetError(MS_GEOSERR,
             "GEOS support enabled, but old version lacks GEOSLength().",
             "msGEOSLength()");
  return -1;
#else
  msSetError(MS_GEOSERR, "GEOS support is not available.", "msGEOSLength()");
  return -1;
#endif
}

double msGEOSDistance(shapeObj *shape1, shapeObj *shape2) {
#ifdef USE_GEOS
  GEOSGeom g1, g2;
  double distance;
  int result;
  GEOSContextHandle_t handle = msGetGeosContextHandle();

  if (!shape1 || !shape2)
    return -1;

  if (!shape1->geometry) /* if no geometry for shape1 then build one */
    shape1->geometry = (GEOSGeom)msGEOSShape2Geometry(shape1);
  g1 = (GEOSGeom)shape1->geometry;
  if (!g1)
    return -1;

  if (!shape2->geometry) /* if no geometry for shape2 then build one */
    shape2->geometry = (GEOSGeom)msGEOSShape2Geometry(shape2);
  g2 = (GEOSGeom)shape2->geometry;
  if (!g2)
    return -1;

  result = GEOSDistance_r(handle, g1, g2, &distance);
  return ((result == 0) ? -1 : distance);
#else
  return msDistanceShapeToShape(
      shape1, shape2); /* fall back on brute force method (for MapScript) */
#endif
}
