/src/postgis/liblwgeom/lwgeom.c

Source
/**********************************************************************
 *
 * PostGIS - Spatial Types for PostgreSQL
 * http://postgis.net
 *
 * PostGIS is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 2 of the License, or
 * (at your option) any later version.
 *
 * PostGIS is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with PostGIS.  If not, see <http://www.gnu.org/licenses/>.
 *
 **********************************************************************
 *
 * Copyright (C) 2001-2006 Refractions Research Inc.
 * Copyright (C) 2017-2018 Daniel Baston <dbaston@gmail.com>
 *
 **********************************************************************/


#include <stdio.h>
#include <stdlib.h>
#include <stdarg.h>
#include <string.h>

#include "liblwgeom_internal.h"
#include "lwgeom_log.h"

#define out_stack_size 32

static char
lwnurbscurve_weight_same(const LWNURBSCURVE *curve1, const LWNURBSCURVE *curve2)
{
  uint32_t npoints = curve1->points->npoints;

  if ((curve1->nweights && curve1->nweights != npoints) ||
      (curve2->nweights && curve2->nweights != npoints))
    return LW_FALSE;

  for (uint32_t i = 0; i < npoints; i++)
  {
    double weight1 = curve1->weights ? curve1->weights[i] : 1.0;
    double weight2 = curve2->weights ? curve2->weights[i] : 1.0;
    if (!FP_EQUALS(weight1, weight2))
      return LW_FALSE;
  }

  return LW_TRUE;
}

static double
lwnurbscurve_uniform_knot_value(uint32_t degree, uint32_t npoints, uint32_t idx)
{
  uint32_t nknots = npoints + degree + 1;
  uint32_t internal_knots = nknots - 2 * (degree + 1);

  if (idx <= degree)
    return 0.0;
  if (idx >= nknots - degree - 1)
    return 1.0;

  return (double)(idx - degree) / (internal_knots + 1);
}

static char
lwnurbscurve_knot_same(const LWNURBSCURVE *curve1, const LWNURBSCURVE *curve2)
{
  uint32_t npoints = curve1->points->npoints;
  uint32_t nknots;

  if (npoints < curve1->degree + 1)
    return LW_FALSE;

  nknots = npoints + curve1->degree + 1;
  if ((curve1->nknots && curve1->nknots != nknots) ||
      (curve2->nknots && curve2->nknots != nknots))
    return LW_FALSE;

  for (uint32_t i = 0; i < nknots; i++)
  {
    double knot1 = curve1->knots ? curve1->knots[i] :
      lwnurbscurve_uniform_knot_value(curve1->degree, npoints, i);
    double knot2 = curve2->knots ? curve2->knots[i] :
      lwnurbscurve_uniform_knot_value(curve2->degree, npoints, i);
    if (!FP_EQUALS(knot1, knot2))
      return LW_FALSE;
  }

  return LW_TRUE;
}

static char
lwnurbscurve_same(const LWNURBSCURVE *curve1, const LWNURBSCURVE *curve2)
{
  if (lwnurbscurve_is_empty(curve1) && lwnurbscurve_is_empty(curve2))
    return LW_TRUE;

  if (curve1->degree != curve2->degree)
    return LW_FALSE;
  if (!curve1->points || !curve2->points || !ptarray_same(curve1->points, curve2->points))
    return LW_FALSE;
  if (!lwnurbscurve_weight_same(curve1, curve2))
    return LW_FALSE;
  if (!lwnurbscurve_knot_same(curve1, curve2))
    return LW_FALSE;

  return LW_TRUE;
}


/**
 * Force an orientation onto geometries made up of rings,
 * so the rings circle in a particular direction, clockwise
 * or counter-clockwise. The orientation is one of
 * LW_CLOCKWISE or LW_COUNTERCLOCKWISE. Non-polygonal
 * geometries are ignored and returned untouched.
 */
static void
lwgeom_force_orientation(LWGEOM *lwgeom, int orientation)
{
  switch (lwgeom->type)
  {
    case POLYGONTYPE:
      lwpoly_force_orientation(lwgeom_as_lwpoly(lwgeom), orientation);
      return;

    case TRIANGLETYPE:
      lwtriangle_force_orientation(lwgeom_as_lwtriangle(lwgeom), orientation);
      return;

    /* Not handle POLYHEDRALSURFACE and TIN
       as they are supposed to be well oriented */
    case MULTIPOLYGONTYPE:
    case COLLECTIONTYPE:
    {
      LWCOLLECTION* coll = lwgeom_as_lwcollection(lwgeom);
      for (uint32_t i = 0; coll && i < coll->ngeoms; i++)
      {
        lwgeom_force_orientation(coll->geoms[i], orientation);
      }
      return;
    }
  }
}

void
lwgeom_force_clockwise(LWGEOM *lwgeom)
{
  return lwgeom_force_orientation(lwgeom, LW_CLOCKWISE);
}

void
lwgeom_force_counterclockwise(LWGEOM *lwgeom)
{
  return lwgeom_force_orientation(lwgeom, LW_COUNTERCLOCKWISE);
}


/**
 * Check clockwise orientation on LWGEOM polygons
 * returns LW_CLOCKWISE, LW_NONE, LW_COUNTERCLOCKWISE
 * Non-polygonal geometries return LW_NONE
 * Geometries must have *all* rings correctly oriented
 * to return a non-none orientation
 */
int
lwgeom_has_orientation(const LWGEOM *lwgeom, int orientation)
{
  switch (lwgeom->type)
  {
    case POLYGONTYPE:
      return lwpoly_has_orientation(lwgeom_as_lwpoly(lwgeom), orientation);

    case TRIANGLETYPE:
      return lwtriangle_has_orientation(lwgeom_as_lwtriangle(lwgeom), orientation);

    case MULTIPOLYGONTYPE:
    case COLLECTIONTYPE:
    {
      LWCOLLECTION* coll = lwgeom_as_lwcollection(lwgeom);
      for (uint32_t i = 0; coll && i < coll->ngeoms; i++)
      {
        if(!lwgeom_has_orientation(coll->geoms[i], orientation))
          return LW_FALSE;
      }
      return LW_TRUE;
    }

    default:
      return LW_TRUE;
  }
  return LW_FALSE;
}

LWGEOM *
lwgeom_reverse(const LWGEOM *geom)
{
  LWGEOM *geomout = lwgeom_clone_deep(geom);
  lwgeom_reverse_in_place(geomout);
  return geomout;
}

/** Reverse vertex order of LWGEOM **/
void
lwgeom_reverse_in_place(LWGEOM *geom)
{
  uint32_t i;
  LWCOLLECTION *col;
  if (!geom)
    return;

  switch (geom->type)
  {
    case MULTIPOINTTYPE:
    case POINTTYPE:
    {
      return;
    }
    case TRIANGLETYPE:
    case CIRCSTRINGTYPE:
    case LINETYPE:
    {
      LWLINE *line = (LWLINE *)(geom);
      ptarray_reverse_in_place(line->points);
      return;
    }
    case POLYGONTYPE:
    {
      LWPOLY *poly = (LWPOLY *)(geom);
      if (!poly->rings)
        return;
      uint32_t r;
      for (r = 0; r < poly->nrings; r++)
        ptarray_reverse_in_place(poly->rings[r]);
      return;
    }
    /* CompoundCurve needs to also reverse the sub-geometries */
    /* so that the end-points remain coincident */
    case COMPOUNDTYPE:
    {
      uint32_t ngeoms;
      col = (LWCOLLECTION *)(geom);
      if (!col->geoms)
        return;
      ngeoms = col->ngeoms;
      for (i=0; i<ngeoms; i++)
        lwgeom_reverse_in_place(col->geoms[i]);
      for (i=0; i<col->ngeoms/2; i++) {
        LWGEOM* tmp = col->geoms[i];
        col->geoms[i] = col->geoms[ngeoms-i-1];
        col->geoms[ngeoms-i-1] = tmp;
      }
      return;
    }
    case MULTICURVETYPE:
    case MULTILINETYPE:
    case MULTIPOLYGONTYPE:
    case MULTISURFACETYPE:
    case POLYHEDRALSURFACETYPE:
    case TINTYPE:
    case COLLECTIONTYPE:
    case CURVEPOLYTYPE:
    {
      col = (LWCOLLECTION *)(geom);
      if (!col->geoms)
        return;
      for (i=0; i<col->ngeoms; i++)
        lwgeom_reverse_in_place(col->geoms[i]);
      return;
    }
    default:
    {
      lwerror("%s: Unknown geometry type: %s", __func__, lwtype_name(geom->type));
      return;
    }

  }
}

LWLINE *
lwgeom_as_lwline(const LWGEOM *lwgeom)
{
  if ( lwgeom == NULL ) return NULL;
  if ( lwgeom->type == LINETYPE )
    return (LWLINE *)lwgeom;
  else return NULL;
}

LWCIRCSTRING *
lwgeom_as_lwcircstring(const LWGEOM *lwgeom)
{
  if ( lwgeom == NULL ) return NULL;
  if ( lwgeom->type == CIRCSTRINGTYPE )
    return (LWCIRCSTRING *)lwgeom;
  else return NULL;
}

LWCOMPOUND *
lwgeom_as_lwcompound(const LWGEOM *lwgeom)
{
  if ( lwgeom == NULL ) return NULL;
  if ( lwgeom->type == COMPOUNDTYPE )
    return (LWCOMPOUND *)lwgeom;
  else return NULL;
}

LWCURVEPOLY *
lwgeom_as_lwcurvepoly(const LWGEOM *lwgeom)
{
  if ( lwgeom == NULL ) return NULL;
  if ( lwgeom->type == CURVEPOLYTYPE )
    return (LWCURVEPOLY *)lwgeom;
  else return NULL;
}

LWPOLY *
lwgeom_as_lwpoly(const LWGEOM *lwgeom)
{
  if ( lwgeom == NULL ) return NULL;
  if ( lwgeom->type == POLYGONTYPE )
    return (LWPOLY *)lwgeom;
  else return NULL;
}

LWTRIANGLE *
lwgeom_as_lwtriangle(const LWGEOM *lwgeom)
{
  if ( lwgeom == NULL ) return NULL;
  if ( lwgeom->type == TRIANGLETYPE )
    return (LWTRIANGLE *)lwgeom;
  else return NULL;
}

LWCOLLECTION *
lwgeom_as_lwcollection(const LWGEOM *lwgeom)
{
  if ( lwgeom == NULL ) return NULL;
  if ( lwgeom_is_collection(lwgeom) )
    return (LWCOLLECTION*)lwgeom;
  else return NULL;
}

LWMPOINT *
lwgeom_as_lwmpoint(const LWGEOM *lwgeom)
{
  if ( lwgeom == NULL ) return NULL;
  if ( lwgeom->type == MULTIPOINTTYPE )
    return (LWMPOINT *)lwgeom;
  else return NULL;
}

LWMLINE *
lwgeom_as_lwmline(const LWGEOM *lwgeom)
{
  if ( lwgeom == NULL ) return NULL;
  if ( lwgeom->type == MULTILINETYPE )
    return (LWMLINE *)lwgeom;
  else return NULL;
}

LWMPOLY *
lwgeom_as_lwmpoly(const LWGEOM *lwgeom)
{
  if ( lwgeom == NULL ) return NULL;
  if ( lwgeom->type == MULTIPOLYGONTYPE )
    return (LWMPOLY *)lwgeom;
  else return NULL;
}

LWPSURFACE *
lwgeom_as_lwpsurface(const LWGEOM *lwgeom)
{
  if ( lwgeom->type == POLYHEDRALSURFACETYPE )
    return (LWPSURFACE *)lwgeom;
  else return NULL;
}

LWTIN *
lwgeom_as_lwtin(const LWGEOM *lwgeom)
{
  if ( lwgeom->type == TINTYPE )
    return (LWTIN *)lwgeom;
  else return NULL;
}

LWGEOM *lwtin_as_lwgeom(const LWTIN *obj)
{
  return (LWGEOM *)obj;
}

LWGEOM *lwpsurface_as_lwgeom(const LWPSURFACE *obj)
{
  return (LWGEOM *)obj;
}

LWGEOM *lwmpoly_as_lwgeom(const LWMPOLY *obj)
{
  if ( obj == NULL ) return NULL;
  return (LWGEOM *)obj;
}
LWGEOM *lwmline_as_lwgeom(const LWMLINE *obj)
{
  if ( obj == NULL ) return NULL;
  return (LWGEOM *)obj;
}
LWGEOM *lwmpoint_as_lwgeom(const LWMPOINT *obj)
{
  if ( obj == NULL ) return NULL;
  return (LWGEOM *)obj;
}
LWGEOM *lwcollection_as_lwgeom(const LWCOLLECTION *obj)
{
  if ( obj == NULL ) return NULL;
  return (LWGEOM *)obj;
}
LWGEOM *lwcircstring_as_lwgeom(const LWCIRCSTRING *obj)
{
  if ( obj == NULL ) return NULL;
  return (LWGEOM *)obj;
}
LWGEOM *lwcurvepoly_as_lwgeom(const LWCURVEPOLY *obj)
{
  if ( obj == NULL ) return NULL;
  return (LWGEOM *)obj;
}
LWGEOM *lwcompound_as_lwgeom(const LWCOMPOUND *obj)
{
  if ( obj == NULL ) return NULL;
  return (LWGEOM *)obj;
}
LWGEOM *lwpoly_as_lwgeom(const LWPOLY *obj)
{
  if ( obj == NULL ) return NULL;
  return (LWGEOM *)obj;
}
LWGEOM *lwtriangle_as_lwgeom(const LWTRIANGLE *obj)
{
  if ( obj == NULL ) return NULL;
  return (LWGEOM *)obj;
}
LWGEOM *lwline_as_lwgeom(const LWLINE *obj)
{
  if ( obj == NULL ) return NULL;
  return (LWGEOM *)obj;
}
LWGEOM *lwpoint_as_lwgeom(const LWPOINT *obj)
{
  if ( obj == NULL ) return NULL;
  return (LWGEOM *)obj;
}


/**
** Look-up for the correct MULTI* type promotion for singleton types.
*/
uint8_t MULTITYPE[NUMTYPES] =
{
  0,
  MULTIPOINTTYPE,        /*  1 */
  MULTILINETYPE,         /*  2 */
  MULTIPOLYGONTYPE,      /*  3 */
  0,0,0,0,
  MULTICURVETYPE,        /*  8 */
  MULTICURVETYPE,        /*  9 */
  MULTISURFACETYPE,      /* 10 */
  POLYHEDRALSURFACETYPE, /* 11 */
  0, 0,
  TINTYPE,               /* 14 */
  0,
  MULTICURVETYPE         /* 16 */
};

uint8_t lwtype_multitype(uint8_t type)
{
  if (type >= NUMTYPES) return 0;
  return MULTITYPE[type];
}

/**
* Create a new LWGEOM of the appropriate MULTI* type.
*/
LWGEOM *
lwgeom_as_multi(const LWGEOM *lwgeom)
{
  LWGEOM **ogeoms;
  LWGEOM *ogeom = NULL;
  GBOX *box = NULL;
  int type;

  type = lwgeom->type;

  /*
  * PolyhedralSurface and MultiPolygon have identical structures
  * (both are collections of LWPOLY). Convert PolyhedralSurface
  * to MultiPolygon by cloning and changing the type.
  */
  if (type == POLYHEDRALSURFACETYPE)
  {
    ogeom = (LWGEOM *)lwcollection_clone((LWCOLLECTION *)lwgeom);
    ogeom->type = MULTIPOLYGONTYPE;
    return ogeom;
  }

  /*
  * TIN to MultiPolygon: convert each triangle to a polygon.
  * LWTRIANGLE and LWLINE have compatible memory layouts,
  * so we can use lwpoly_from_lwlines to create polygons.
  */
  if (type == TINTYPE)
  {
    uint32_t i;
    LWCOLLECTION *col = (LWCOLLECTION *)lwgeom;
    LWMPOLY *mpoly = lwmpoly_construct_empty(lwgeom->srid,
      FLAGS_GET_Z(lwgeom->flags),
      FLAGS_GET_M(lwgeom->flags));

    for (i = 0; i < col->ngeoms; i++)
    {
      LWPOLY *poly = lwpoly_from_lwlines((LWLINE *)col->geoms[i], 0, NULL);
      lwmpoly_add_lwpoly(mpoly, poly);
    }

    if (lwgeom->bbox)
      mpoly->bbox = gbox_clone(lwgeom->bbox);

    return lwmpoly_as_lwgeom(mpoly);
  }

  if ( ! MULTITYPE[type] ) return lwgeom_clone(lwgeom);

  if( lwgeom_is_empty(lwgeom) )
  {
    ogeom = (LWGEOM *)lwcollection_construct_empty(
      MULTITYPE[type],
      lwgeom->srid,
      FLAGS_GET_Z(lwgeom->flags),
      FLAGS_GET_M(lwgeom->flags)
    );
  }
  else
  {
    ogeoms = lwalloc(sizeof(LWGEOM*));
    ogeoms[0] = lwgeom_clone(lwgeom);

    /* Sub-geometries are not allowed to have bboxes or SRIDs, move the bbox to the collection */
    box = ogeoms[0]->bbox;
    ogeoms[0]->bbox = NULL;
    ogeoms[0]->srid = SRID_UNKNOWN;

    ogeom = (LWGEOM *)lwcollection_construct(MULTITYPE[type], lwgeom->srid, box, 1, ogeoms);
  }

  return ogeom;
}

/**
* Create a new LWGEOM of the appropriate CURVE* type.
*/
LWGEOM *
lwgeom_as_curve(const LWGEOM *lwgeom)
{
  LWGEOM *ogeom;
  int type = lwgeom->type;
  /*
  int hasz = FLAGS_GET_Z(lwgeom->flags);
  int hasm = FLAGS_GET_M(lwgeom->flags);
  int32_t srid = lwgeom->srid;
  */

  switch(type)
  {
    case LINETYPE:
      /* turn to COMPOUNDCURVE */
      ogeom = (LWGEOM*)lwcompound_construct_from_lwline((LWLINE*)lwgeom);
      break;
    case POLYGONTYPE:
      ogeom = (LWGEOM*)lwcurvepoly_construct_from_lwpoly(lwgeom_as_lwpoly(lwgeom));
      break;
    case MULTILINETYPE:
      /* turn to MULTICURVE */
      ogeom = lwgeom_clone(lwgeom);
      ogeom->type = MULTICURVETYPE;
      break;
    case MULTIPOLYGONTYPE:
      /* turn to MULTISURFACE */
      ogeom = lwgeom_clone(lwgeom);
      ogeom->type = MULTISURFACETYPE;
      break;
    case COLLECTIONTYPE:
    default:
      ogeom = lwgeom_clone(lwgeom);
      break;
  }

  /* TODO: copy bbox from input geom ? */

  return ogeom;
}


/**
* Free the containing LWGEOM and the associated BOX. Leave the underlying
* geoms/points/point objects intact. Useful for functions that are stripping
* out subcomponents of complex objects, or building up new temporary objects
* on top of subcomponents.
*/
void
lwgeom_release(LWGEOM *lwgeom)
{
  if ( ! lwgeom )
    lwerror("lwgeom_release: someone called on 0x0");

  LWDEBUGF(3, "releasing type %s", lwtype_name(lwgeom->type));

  /* Drop bounding box (always a copy) */
  if ( lwgeom->bbox )
  {
    LWDEBUGF(3, "lwgeom_release: releasing bbox. %p", lwgeom->bbox);
    lwfree(lwgeom->bbox);
  }
  lwfree(lwgeom);

}


/* @brief Clone LWGEOM object. Serialized point lists are not copied.
 *
 * @see ptarray_clone
 */
LWGEOM *
lwgeom_clone(const LWGEOM *lwgeom)
{
  LWDEBUGF(2, "lwgeom_clone called with %p, %s",
           lwgeom, lwtype_name(lwgeom->type));

  switch (lwgeom->type)
  {
  case POINTTYPE:
    return (LWGEOM *)lwpoint_clone((LWPOINT *)lwgeom);
  case LINETYPE:
    return (LWGEOM *)lwline_clone((LWLINE *)lwgeom);
  case CIRCSTRINGTYPE:
    return (LWGEOM *)lwcircstring_clone((LWCIRCSTRING *)lwgeom);
  case NURBSCURVETYPE:
    return (LWGEOM *)lwnurbscurve_clone_deep((LWNURBSCURVE *)lwgeom);
  case POLYGONTYPE:
    return (LWGEOM *)lwpoly_clone((LWPOLY *)lwgeom);
  case TRIANGLETYPE:
    return (LWGEOM *)lwtriangle_clone((LWTRIANGLE *)lwgeom);
  case COMPOUNDTYPE:
  case CURVEPOLYTYPE:
  case MULTICURVETYPE:
  case MULTISURFACETYPE:
  case MULTIPOINTTYPE:
  case MULTILINETYPE:
  case MULTIPOLYGONTYPE:
  case POLYHEDRALSURFACETYPE:
  case TINTYPE:
  case COLLECTIONTYPE:
    return (LWGEOM *)lwcollection_clone((LWCOLLECTION *)lwgeom);
  default:
    lwerror("lwgeom_clone: Unknown geometry type: %s", lwtype_name(lwgeom->type));
    return NULL;
  }
}

/**
 * Create a deep copy of an LWGEOM.
 *
 * Performs a full (deep) clone of the provided geometry, including copying
 * internal vertex/point arrays so the returned geometry is independent of the
 * original. Dispatches to type-specific deep-clone helpers (e.g. line, polygon,
 * collection, NURBS curve).
 *
 * @param lwgeom Source geometry to clone. Must be non-NULL.
 * @returns Pointer to a newly allocated LWGEOM that is a deep clone of
 *          `lwgeom`, or NULL if `lwgeom->type` is unrecognized (an error is
 *          reported in that case).
 */
LWGEOM *
lwgeom_clone_deep(const LWGEOM *lwgeom)
{
  LWDEBUGF(2, "lwgeom_clone called with %p, %s",
           lwgeom, lwtype_name(lwgeom->type));

  switch (lwgeom->type)
  {
  case POINTTYPE:
  case LINETYPE:
  case CIRCSTRINGTYPE:
  case TRIANGLETYPE:
    return (LWGEOM *)lwline_clone_deep((LWLINE *)lwgeom);
  case POLYGONTYPE:
    return (LWGEOM *)lwpoly_clone_deep((LWPOLY *)lwgeom);
  case COMPOUNDTYPE:
  case CURVEPOLYTYPE:
  case MULTICURVETYPE:
  case MULTISURFACETYPE:
  case MULTIPOINTTYPE:
  case MULTILINETYPE:
  case MULTIPOLYGONTYPE:
  case POLYHEDRALSURFACETYPE:
  case TINTYPE:
  case COLLECTIONTYPE:
    return (LWGEOM *)lwcollection_clone_deep((LWCOLLECTION *)lwgeom);
  case NURBSCURVETYPE:
    return (LWGEOM *)lwnurbscurve_clone_deep((LWNURBSCURVE *)lwgeom);
  default:
    lwerror("lwgeom_clone_deep: Unknown geometry type: %s", lwtype_name(lwgeom->type));
    return NULL;
  }
}


/**
 * Return an allocated string
 */
char*
lwgeom_to_ewkt(const LWGEOM *lwgeom)
{
  char* wkt = NULL;
  size_t wkt_size = 0;

  wkt = lwgeom_to_wkt(lwgeom, WKT_EXTENDED, 12, &wkt_size);

  if ( ! wkt )
  {
    lwerror("Error writing geom %p to WKT", (void *)lwgeom);
  }

  return wkt;
}

/**
 * @brief geom1 same as geom2
 *    iff
 *        + have same type
 *      + have same # objects
 *        + have same bvol
 *        + each object in geom1 has a corresponding object in geom2 (see above)
 *  @param lwgeom1
 *  @param lwgeom2
 */
char
lwgeom_same(const LWGEOM *lwgeom1, const LWGEOM *lwgeom2)
{
  LWDEBUGF(2, "lwgeom_same(%s, %s) called",
           lwtype_name(lwgeom1->type),
           lwtype_name(lwgeom2->type));

  if ( lwgeom1->type != lwgeom2->type )
  {
    LWDEBUG(3, " type differ");

    return LW_FALSE;
  }

  if ( FLAGS_GET_ZM(lwgeom1->flags) != FLAGS_GET_ZM(lwgeom2->flags) )
  {
    LWDEBUG(3, " ZM flags differ");

    return LW_FALSE;
  }

  /* Check boxes if both already computed  */
  if ( lwgeom1->bbox && lwgeom2->bbox )
  {
    /*lwnotice("bbox1:%p, bbox2:%p", lwgeom1->bbox, lwgeom2->bbox);*/
    if ( ! gbox_same(lwgeom1->bbox, lwgeom2->bbox) )
    {
      LWDEBUG(3, " bounding boxes differ");

      return LW_FALSE;
    }
  }

  /* geoms have same type, invoke type-specific function */
  switch (lwgeom1->type)
  {
  case POINTTYPE:
    return lwpoint_same((LWPOINT *)lwgeom1,
                        (LWPOINT *)lwgeom2);
  case LINETYPE:
    return lwline_same((LWLINE *)lwgeom1,
                       (LWLINE *)lwgeom2);
  case POLYGONTYPE:
    return lwpoly_same((LWPOLY *)lwgeom1,
                       (LWPOLY *)lwgeom2);
  case TRIANGLETYPE:
    return lwtriangle_same((LWTRIANGLE *)lwgeom1,
                           (LWTRIANGLE *)lwgeom2);
  case CIRCSTRINGTYPE:
    return lwcircstring_same((LWCIRCSTRING *)lwgeom1,
           (LWCIRCSTRING *)lwgeom2);
  case NURBSCURVETYPE:
    return lwnurbscurve_same((LWNURBSCURVE *)lwgeom1,
                             (LWNURBSCURVE *)lwgeom2);
  case MULTIPOINTTYPE:
  case MULTILINETYPE:
  case MULTIPOLYGONTYPE:
  case MULTICURVETYPE:
  case MULTISURFACETYPE:
  case COMPOUNDTYPE:
  case CURVEPOLYTYPE:
  case POLYHEDRALSURFACETYPE:
  case TINTYPE:
  case COLLECTIONTYPE:
    return lwcollection_same((LWCOLLECTION *)lwgeom1,
                             (LWCOLLECTION *)lwgeom2);
  default:
    lwerror("lwgeom_same: unsupported geometry type: %s",
            lwtype_name(lwgeom1->type));
    return LW_FALSE;
  }

}

int
lwpoint_inside_circle(const LWPOINT *p, double cx, double cy, double rad)
{
  const POINT2D *pt;
  POINT2D center;

  if ( ! p || ! p->point )
    return LW_FALSE;

  pt = getPoint2d_cp(p->point, 0);

  center.x = cx;
  center.y = cy;

  if ( distance2d_pt_pt(pt, &center) < rad )
    return LW_TRUE;

  return LW_FALSE;
}

void
lwgeom_drop_bbox(LWGEOM *lwgeom)
{
  if ( lwgeom->bbox ) lwfree(lwgeom->bbox);
  lwgeom->bbox = NULL;
  FLAGS_SET_BBOX(lwgeom->flags, 0);
}

/**
 * Ensure there's a box in the LWGEOM.
 * If the box is already there just return,
 * else compute it.
 */
void
lwgeom_add_bbox(LWGEOM *lwgeom)
{
  /* an empty LWGEOM has no bbox */
  if ( lwgeom_is_empty(lwgeom) ) return;

  if ( lwgeom->bbox ) return;
  FLAGS_SET_BBOX(lwgeom->flags, 1);
  lwgeom->bbox = gbox_new(lwgeom->flags);
  lwgeom_calculate_gbox(lwgeom, lwgeom->bbox);
}

void
lwgeom_refresh_bbox(LWGEOM *lwgeom)
{
  lwgeom_drop_bbox(lwgeom);
  lwgeom_add_bbox(lwgeom);
}

void
lwgeom_add_bbox_deep(LWGEOM *lwgeom, GBOX *gbox)
{
  if ( lwgeom_is_empty(lwgeom) ) return;

  FLAGS_SET_BBOX(lwgeom->flags, 1);

  if ( ! ( gbox || lwgeom->bbox ) )
  {
    lwgeom->bbox = gbox_new(lwgeom->flags);
    lwgeom_calculate_gbox(lwgeom, lwgeom->bbox);
  }
  else if ( gbox && ! lwgeom->bbox )
  {
    lwgeom->bbox = gbox_clone(gbox);
  }

  if ( lwgeom_is_collection(lwgeom) )
  {
    uint32_t i;
    LWCOLLECTION *lwcol = (LWCOLLECTION*)lwgeom;

    for ( i = 0; i < lwcol->ngeoms; i++ )
    {
      lwgeom_add_bbox_deep(lwcol->geoms[i], lwgeom->bbox);
    }
  }
}

const GBOX *
lwgeom_get_bbox(const LWGEOM *lwg)
{
  /* add it if not already there */
  lwgeom_add_bbox((LWGEOM *)lwg);
  return lwg->bbox;
}


/**
* Calculate the gbox for this geometry, a cartesian box or
* geodetic box, depending on how it is flagged.
*/
int lwgeom_calculate_gbox(const LWGEOM *lwgeom, GBOX *gbox)
{
  gbox->flags = lwgeom->flags;
  if( FLAGS_GET_GEODETIC(lwgeom->flags) )
    return lwgeom_calculate_gbox_geodetic(lwgeom, gbox);
  else
    return lwgeom_calculate_gbox_cartesian(lwgeom, gbox);
}

void
lwgeom_drop_srid(LWGEOM *lwgeom)
{
  lwgeom->srid = SRID_UNKNOWN; /* TODO: To be changed to SRID_UNKNOWN */
}

LWGEOM *
lwgeom_segmentize2d(const LWGEOM *lwgeom, double dist)
{
  switch (lwgeom->type)
  {
  case LINETYPE:
    return (LWGEOM *)lwline_segmentize2d((LWLINE *)lwgeom,
                                         dist);
  case POLYGONTYPE:
    return (LWGEOM *)lwpoly_segmentize2d((LWPOLY *)lwgeom,
                                         dist);
  case MULTILINETYPE:
  case MULTIPOLYGONTYPE:
  case COLLECTIONTYPE:
    return (LWGEOM *)lwcollection_segmentize2d(
               (LWCOLLECTION *)lwgeom, dist);

  default:
    return lwgeom_clone(lwgeom);
  }
}

LWGEOM*
lwgeom_force_2d(const LWGEOM *geom)
{
  return lwgeom_force_dims(geom, 0, 0, 0, 0);
}

LWGEOM*
lwgeom_force_3dz(const LWGEOM *geom, double zval)
{
  return lwgeom_force_dims(geom, 1, 0, zval, 0);
}

LWGEOM*
lwgeom_force_3dm(const LWGEOM *geom, double mval)
{
  return lwgeom_force_dims(geom, 0, 1, 0, mval);
}

LWGEOM*
lwgeom_force_4d(const LWGEOM *geom, double zval, double mval)
{
  return lwgeom_force_dims(geom, 1, 1, zval, mval);
}

LWGEOM*
lwgeom_force_dims(const LWGEOM *geom, int hasz, int hasm, double zval, double mval)
{
  if (!geom)
    return NULL;
  switch(geom->type)
  {
    case POINTTYPE:
      return lwpoint_as_lwgeom(lwpoint_force_dims((LWPOINT*)geom, hasz, hasm, zval, mval));
    case CIRCSTRINGTYPE:
    case LINETYPE:
    case TRIANGLETYPE:
      return lwline_as_lwgeom(lwline_force_dims((LWLINE*)geom, hasz, hasm, zval, mval));
    case POLYGONTYPE:
      return lwpoly_as_lwgeom(lwpoly_force_dims((LWPOLY*)geom, hasz, hasm, zval, mval));
    case COMPOUNDTYPE:
    case CURVEPOLYTYPE:
    case MULTICURVETYPE:
    case MULTISURFACETYPE:
    case MULTIPOINTTYPE:
    case MULTILINETYPE:
    case MULTIPOLYGONTYPE:
    case POLYHEDRALSURFACETYPE:
    case TINTYPE:
    case COLLECTIONTYPE:
      return lwcollection_as_lwgeom(lwcollection_force_dims((LWCOLLECTION*)geom, hasz, hasm, zval, mval));
    case NURBSCURVETYPE:
      return lwnurbscurve_as_lwgeom(lwnurbscurve_force_dims((LWNURBSCURVE*)geom, hasz, hasm, zval, mval));
    default:
      lwerror("lwgeom_force_2d: unsupported geom type: %s", lwtype_name(geom->type));
      return NULL;
  }
}

LWGEOM*
lwgeom_force_sfs(LWGEOM *geom, int version)
{
  LWCOLLECTION *col;
  uint32_t i;
  LWGEOM *g;

  /* SFS 1.2 version */
  if (version == 120)
  {
    switch(geom->type)
    {
      /* SQL/MM types */
      case CIRCSTRINGTYPE:
      case COMPOUNDTYPE:
      case CURVEPOLYTYPE:
      case MULTICURVETYPE:
      case MULTISURFACETYPE:
        return lwgeom_stroke(geom, 32);

      case COLLECTIONTYPE:
        col = (LWCOLLECTION*)geom;
        for ( i = 0; i < col->ngeoms; i++ )
          col->geoms[i] = lwgeom_force_sfs((LWGEOM*)col->geoms[i], version);

        return lwcollection_as_lwgeom((LWCOLLECTION*)geom);

      default:
        return (LWGEOM *)geom;
    }
  }


  /* SFS 1.1 version */
  switch(geom->type)
  {
    /* SQL/MM types */
    case CIRCSTRINGTYPE:
    case COMPOUNDTYPE:
    case CURVEPOLYTYPE:
    case MULTICURVETYPE:
    case MULTISURFACETYPE:
      return lwgeom_stroke(geom, 32);

    /* SFS 1.2 types */
    case TRIANGLETYPE:
      g = lwpoly_as_lwgeom(lwpoly_from_lwlines((LWLINE*)geom, 0, NULL));
      lwgeom_free(geom);
      return g;

    case TINTYPE:
      col = (LWCOLLECTION*) geom;
      for ( i = 0; i < col->ngeoms; i++ )
      {
        g = lwpoly_as_lwgeom(lwpoly_from_lwlines((LWLINE*)col->geoms[i], 0, NULL));
        lwgeom_free(col->geoms[i]);
        col->geoms[i] = g;
      }
      col->type = COLLECTIONTYPE;
      return lwmpoly_as_lwgeom((LWMPOLY*)geom);

    case POLYHEDRALSURFACETYPE:
      geom->type = COLLECTIONTYPE;
      return (LWGEOM *)geom;

    /* Collection */
    case COLLECTIONTYPE:
      col = (LWCOLLECTION*)geom;
      for ( i = 0; i < col->ngeoms; i++ )
        col->geoms[i] = lwgeom_force_sfs((LWGEOM*)col->geoms[i], version);

      return lwcollection_as_lwgeom((LWCOLLECTION*)geom);

    default:
      return (LWGEOM *)geom;
  }
}

int32_t
lwgeom_get_srid(const LWGEOM *geom)
{
  if ( ! geom ) return SRID_UNKNOWN;
  return geom->srid;
}

int
lwgeom_has_z(const LWGEOM *geom)
{
  if ( ! geom ) return LW_FALSE;
  return FLAGS_GET_Z(geom->flags);
}

int
lwgeom_has_m(const LWGEOM *geom)
{
  if ( ! geom ) return LW_FALSE;
  return FLAGS_GET_M(geom->flags);
}

int
lwgeom_is_solid(const LWGEOM *geom)
{
  if ( ! geom ) return LW_FALSE;
  return FLAGS_GET_SOLID(geom->flags);
}

int
lwgeom_ndims(const LWGEOM *geom)
{
  if ( ! geom ) return 0;
  return FLAGS_NDIMS(geom->flags);
}



void
lwgeom_set_geodetic(LWGEOM *geom, int value)
{
  LWPOINT *pt;
  LWLINE *ln;
  LWPOLY *ply;
  LWCOLLECTION *col;
  uint32_t i;

  FLAGS_SET_GEODETIC(geom->flags, value);
  if ( geom->bbox )
    FLAGS_SET_GEODETIC(geom->bbox->flags, value);

  switch(geom->type)
  {
    case POINTTYPE:
      pt = (LWPOINT*)geom;
      if ( pt->point )
        FLAGS_SET_GEODETIC(pt->point->flags, value);
      break;
    case LINETYPE:
      ln = (LWLINE*)geom;
      if ( ln->points )
        FLAGS_SET_GEODETIC(ln->points->flags, value);
      break;
    case POLYGONTYPE:
      ply = (LWPOLY*)geom;
      for ( i = 0; i < ply->nrings; i++ )
        FLAGS_SET_GEODETIC(ply->rings[i]->flags, value);
      break;
    case MULTIPOINTTYPE:
    case MULTILINETYPE:
    case MULTIPOLYGONTYPE:
    case COLLECTIONTYPE:
      col = (LWCOLLECTION*)geom;
      for ( i = 0; i < col->ngeoms; i++ )
        lwgeom_set_geodetic(col->geoms[i], value);
      break;
    default:
      lwerror("lwgeom_set_geodetic: unsupported geom type: %s", lwtype_name(geom->type));
      return;
  }
}

void
lwgeom_longitude_shift(LWGEOM *lwgeom)
{
  uint32_t i;
  switch (lwgeom->type)
  {
    LWPOINT *point;
    LWLINE *line;
    LWPOLY *poly;
    LWTRIANGLE *triangle;
    LWCOLLECTION *coll;

  case POINTTYPE:
    point = (LWPOINT *)lwgeom;
    ptarray_longitude_shift(point->point);
    return;
  case LINETYPE:
    line = (LWLINE *)lwgeom;
    ptarray_longitude_shift(line->points);
    return;
  case POLYGONTYPE:
    poly = (LWPOLY *)lwgeom;
    for (i=0; i<poly->nrings; i++)
      ptarray_longitude_shift(poly->rings[i]);
    return;
  case TRIANGLETYPE:
    triangle = (LWTRIANGLE *)lwgeom;
    ptarray_longitude_shift(triangle->points);
    return;
  case MULTIPOINTTYPE:
  case MULTILINETYPE:
  case MULTIPOLYGONTYPE:
  case POLYHEDRALSURFACETYPE:
  case TINTYPE:
  case COLLECTIONTYPE:
    coll = (LWCOLLECTION *)lwgeom;
    for (i=0; i<coll->ngeoms; i++)
      lwgeom_longitude_shift(coll->geoms[i]);
    return;
  default:
    lwerror("lwgeom_longitude_shift: unsupported geom type: %s",
            lwtype_name(lwgeom->type));
  }
}

int
lwgeom_is_closed(const LWGEOM *geom)
{
  int type = geom->type;

  if( lwgeom_is_empty(geom) )
    return LW_FALSE;

  /* Test linear types for closure */
  switch (type)
  {
  case LINETYPE:
    return lwline_is_closed((LWLINE*)geom);
  case POLYGONTYPE:
    return lwpoly_is_closed((LWPOLY*)geom);
  case CIRCSTRINGTYPE:
    return lwcircstring_is_closed((LWCIRCSTRING*)geom);
  case NURBSCURVETYPE:
  {
    POINT4D start, end;
    LWPOINT *startpt = lwnurbscurve_evaluate((const LWNURBSCURVE*)geom, 0.0);
    LWPOINT *endpt = lwnurbscurve_evaluate((const LWNURBSCURVE*)geom, 1.0);
    int closed = LW_FALSE;

    if (startpt && endpt &&
        lwpoint_getPoint4d_p(startpt, &start) &&
        lwpoint_getPoint4d_p(endpt, &end))
    {
      if (FLAGS_GET_Z(geom->flags))
        closed = FP_EQUALS(start.x, end.x) && FP_EQUALS(start.y, end.y) && FP_EQUALS(start.z, end.z);
      else
        closed = FP_EQUALS(start.x, end.x) && FP_EQUALS(start.y, end.y);
    }

    lwpoint_free(startpt);
    lwpoint_free(endpt);
    return closed;
  }
  case COMPOUNDTYPE:
    return lwcompound_is_closed((LWCOMPOUND*)geom);
  case TINTYPE:
    return lwtin_is_closed((LWTIN*)geom);
  case POLYHEDRALSURFACETYPE:
    return lwpsurface_is_closed((LWPSURFACE*)geom);
  }

  /* Recurse into collections and see if anything is not closed */
  if ( lwgeom_is_collection(geom) )
  {
    LWCOLLECTION *col = lwgeom_as_lwcollection(geom);
    uint32_t i;
    int closed;
    for ( i = 0; i < col->ngeoms; i++ )
    {
      closed = lwgeom_is_closed(col->geoms[i]);
      if ( ! closed )
        return LW_FALSE;
    }
    return LW_TRUE;
  }

  /* All non-linear non-collection types we will call closed */
  return LW_TRUE;
}

int
lwgeom_is_collection(const LWGEOM *geom)
{
  if( ! geom ) return LW_FALSE;
  return lwtype_is_collection(geom->type);
}

int
lwtype_is_unitary(uint32_t lwtype)
{
  switch (lwtype)
  {
  case POINTTYPE:
  case LINETYPE:
  case POLYGONTYPE:
  case CURVEPOLYTYPE:
  case COMPOUNDTYPE:
  case CIRCSTRINGTYPE:
  case NURBSCURVETYPE:
  case TRIANGLETYPE:
  case POLYHEDRALSURFACETYPE:
  case TINTYPE:
    return LW_TRUE;
    break;

  default:
    return LW_FALSE;
  }
}

int
lwgeom_has_patches(const LWGEOM *geom)
{
  switch (geom->type)
  {
  case TINTYPE:
  case POLYHEDRALSURFACETYPE:
    return LW_TRUE;
    break;
  default:
    return LW_FALSE;
  }
}

int
lwgeom_is_unitary(const LWGEOM *geom)
{
  return lwtype_is_unitary(geom->type);
}

int
lwgeom_has_rings(const LWGEOM *geom)
{
  switch (geom->type)
  {
  case POLYGONTYPE:
  case CURVEPOLYTYPE:
  case TRIANGLETYPE:
    return LW_TRUE;
    break;

  default:
    return LW_FALSE;
  }
}

/**
 * Return TRUE if the geometry is structured as a wrapper on a geoms/ngeoms
 * list of sub-geometries.
 * Use lwgeom_is_unitary to determine if the numgeometries/geometryn accessor
 * pattery makes sense
 */
int
lwtype_is_collection(uint8_t type)
{
  switch (type)
  {
  case MULTIPOINTTYPE:
  case MULTILINETYPE:
  case MULTIPOLYGONTYPE:
  case COLLECTIONTYPE:
  case CURVEPOLYTYPE:
  case COMPOUNDTYPE:
  case MULTICURVETYPE:
  case MULTISURFACETYPE:
  case POLYHEDRALSURFACETYPE:
  case TINTYPE:
    return LW_TRUE;
    break;

  default:
    return LW_FALSE;
  }
}

/**
* Given an lwtype number, what homogeneous collection can hold it?
*/
uint32_t
lwtype_get_collectiontype(uint8_t type)
{
  switch (type)
  {
    case POINTTYPE:
      return MULTIPOINTTYPE;
    case LINETYPE:
      return MULTILINETYPE;
    case POLYGONTYPE:
      return MULTIPOLYGONTYPE;
    case CIRCSTRINGTYPE:
      return MULTICURVETYPE;
    case COMPOUNDTYPE:
      return MULTICURVETYPE;
    case NURBSCURVETYPE:
      return MULTICURVETYPE;
    case CURVEPOLYTYPE:
      return MULTISURFACETYPE;
    case TRIANGLETYPE:
      return TINTYPE;
    default:
      return COLLECTIONTYPE;
  }
}


/**
 * Free an LWGEOM and its associated resources.
 *
 * Dispatches to the appropriate type-specific free function based on lwgeom->type.
 * If `lwgeom` is NULL the function is a no-op. Logs an error if the geometry type
 * is unknown.
 *
 * @param lwgeom Pointer to the LWGEOM to free (may be NULL). After return the
 *               memory for the provided geometry will have been released.
 */
void lwgeom_free(LWGEOM *lwgeom)
{

  /* There's nothing here to free... */
  if( ! lwgeom ) return;

  LWDEBUGF(5,"freeing a %s",lwtype_name(lwgeom->type));

  switch (lwgeom->type)
  {
  case POINTTYPE:
    lwpoint_free((LWPOINT *)lwgeom);
    break;
  case LINETYPE:
    lwline_free((LWLINE *)lwgeom);
    break;
  case POLYGONTYPE:
    lwpoly_free((LWPOLY *)lwgeom);
    break;
  case CIRCSTRINGTYPE:
    lwcircstring_free((LWCIRCSTRING *)lwgeom);
    break;
  case TRIANGLETYPE:
    lwtriangle_free((LWTRIANGLE *)lwgeom);
    break;
  case MULTIPOINTTYPE:
    lwmpoint_free((LWMPOINT *)lwgeom);
    break;
  case MULTILINETYPE:
    lwmline_free((LWMLINE *)lwgeom);
    break;
  case MULTIPOLYGONTYPE:
    lwmpoly_free((LWMPOLY *)lwgeom);
    break;
  case POLYHEDRALSURFACETYPE:
    lwpsurface_free((LWPSURFACE *)lwgeom);
    break;
  case TINTYPE:
    lwtin_free((LWTIN *)lwgeom);
    break;
  case CURVEPOLYTYPE:
  case COMPOUNDTYPE:
  case MULTICURVETYPE:
  case MULTISURFACETYPE:
  case COLLECTIONTYPE:
    lwcollection_free((LWCOLLECTION *)lwgeom);
    break;
  case NURBSCURVETYPE:
    lwnurbscurve_free((LWNURBSCURVE *)lwgeom);
    break;
  default:
    lwerror("lwgeom_free called with unknown type (%d) %s", lwgeom->type, lwtype_name(lwgeom->type));
  }
  return;
}

/**
 * Determine whether a geometry requires a bounding box.
 *
 * Returns LW_TRUE if callers should attach/maintain a bounding box for the
 * provided geometry type and state, or LW_FALSE when a bbox is unnecessary
 * (e.g., point-like or trivially short line geometries).
 *
 * Detailed behavior:
 * - POINTTYPE: returns LW_FALSE.
 * - LINETYPE: returns LW_FALSE for lines with 0–2 vertices, otherwise LW_TRUE.
 * - NURBSCURVETYPE: currently treated as not requiring a bbox (LW_FALSE).
 * - MULTIPOINTTYPE: returns LW_FALSE for a single-member collection, otherwise LW_TRUE.
 * - MULTILINETYPE: returns LW_FALSE when the collection has one member and that member has ≤2 vertices, otherwise LW_TRUE.
 * - All other types: returns LW_TRUE.
 *
 * @param geom Geometry to inspect (must not be NULL).
 * @return LW_TRUE if a bbox is needed, LW_FALSE otherwise.
 */
int lwgeom_needs_bbox(const LWGEOM *geom)
{
  assert(geom);
  if ( geom->type == POINTTYPE )
  {
    return LW_FALSE;
  }
  else if ( geom->type == LINETYPE )
  {
    if ( lwgeom_count_vertices(geom) <= 2 )
      return LW_FALSE;
    else
      return LW_TRUE;
  }
  else if ( geom->type == NURBSCURVETYPE )
  {
    const LWNURBSCURVE *nurbs = (LWNURBSCURVE*)geom;
    /* Empty curves don't need bbox */
    if ( !nurbs->points || nurbs->points->npoints == 0 )
      return LW_FALSE;
    /* Single control point doesn't need bbox */
    else if ( nurbs->points->npoints == 1 )
      return LW_FALSE;
    else
      return LW_TRUE;
  }
  else if ( geom->type == MULTIPOINTTYPE )
  {
    if ( ((LWCOLLECTION*)geom)->ngeoms == 1 )
      return LW_FALSE;
    else
      return LW_TRUE;
  }
  else if ( geom->type == MULTILINETYPE )
  {
    if ( ((LWCOLLECTION*)geom)->ngeoms == 1 && lwgeom_count_vertices(geom) <= 2 )
      return LW_FALSE;
    else
      return LW_TRUE;
  }
  else
  {
    return LW_TRUE;
  }
}

/**
* Count points in an #LWGEOM.
* TODO: Make sure the internal functions don't overflow
*/
uint32_t lwgeom_count_vertices(const LWGEOM *geom)
{
  int result = 0;

  /* Null? Zero. */
  if( ! geom ) return 0;

  LWDEBUGF(4, "lwgeom_count_vertices got type %s",
           lwtype_name(geom->type));

  /* Empty? Zero. */
  if( lwgeom_is_empty(geom) ) return 0;

  switch (geom->type)
  {
  case POINTTYPE:
    result = 1;
    break;
  case TRIANGLETYPE:
  case CIRCSTRINGTYPE:
  case LINETYPE:
    result = lwline_count_vertices((const LWLINE *)geom);
    break;
  case NURBSCURVETYPE:
    {
      const LWNURBSCURVE *nurbs = (const LWNURBSCURVE *)geom;
      result = (nurbs->points != NULL) ? nurbs->points->npoints : 0;
    }
    break;
  case POLYGONTYPE:
    result = lwpoly_count_vertices((const LWPOLY *)geom);
    break;
  case COMPOUNDTYPE:
  case CURVEPOLYTYPE:
  case MULTICURVETYPE:
  case MULTISURFACETYPE:
  case MULTIPOINTTYPE:
  case MULTILINETYPE:
  case MULTIPOLYGONTYPE:
  case POLYHEDRALSURFACETYPE:
  case TINTYPE:
  case COLLECTIONTYPE:
    result = lwcollection_count_vertices((const LWCOLLECTION *)geom);
    break;
  default:
    lwerror("%s: unsupported input geometry type: %s",
            __func__, lwtype_name(geom->type));
    break;
  }
  LWDEBUGF(3, "counted %d vertices", result);
  return result;
}

/**
 * Return the topological dimension of a geometry.
 *
 * Determines the topological dimension for the given LWGEOM:
 * - 0 for point-like geometries,
 * - 1 for curve-like geometries (lines, circular strings, compound/curve collections, NURBS),
 * - 2 for surface-like geometries (polygons, triangles, surfaces, TIN),
 * - 3 for a closed polyhedral surface (volume).
 * For a COLLECTION, returns the maximum dimension among its members.
 *
 * @param geom Pointer to the geometry; may be NULL.
 * @return 0..3 for the geometry's topological dimension, the maximum member dimension for collections,
 *         -1 if geom is NULL or its type is unsupported.
 */
int lwgeom_dimension(const LWGEOM *geom)
{

  /* Null? Zero. */
  if( ! geom ) return -1;

  LWDEBUGF(4, "lwgeom_dimension got type %s",
           lwtype_name(geom->type));

  /* Empty? Zero. */
  /* if( lwgeom_is_empty(geom) ) return 0; */

  switch (geom->type)
  {
  case POINTTYPE:
  case MULTIPOINTTYPE:
    return 0;
  case CIRCSTRINGTYPE:
  case LINETYPE:
  case COMPOUNDTYPE:
  case MULTICURVETYPE:
  case MULTILINETYPE:
  case NURBSCURVETYPE:
    return 1;
  case TRIANGLETYPE:
  case POLYGONTYPE:
  case CURVEPOLYTYPE:
  case MULTISURFACETYPE:
  case MULTIPOLYGONTYPE:
    return 2;
  case TINTYPE:
    /* A 3D TIN surface is dimension 3 only when explicitly marked solid. */
    return (FLAGS_GET_SOLID(geom->flags) && FLAGS_GET_Z(geom->flags)) ? 3 : 2;
  case POLYHEDRALSURFACETYPE:
  {
    /* A closed polyhedral surface contains a volume. */
    int closed = lwpsurface_is_closed((LWPSURFACE*)geom);
    return ( closed ? 3 : 2 );
  }
  case COLLECTIONTYPE:
  {
    int maxdim = 0;
    uint32_t i;
    LWCOLLECTION *col = (LWCOLLECTION*)geom;
    for( i = 0; i < col->ngeoms; i++ )
    {
      int dim = lwgeom_dimension(col->geoms[i]);
      maxdim = ( dim > maxdim ? dim : maxdim );
    }
    return maxdim;
  }
  default:
    lwerror("%s: unsupported input geometry type: %s",
            __func__, lwtype_name(geom->type));
  }
  return -1;
}

/**
* Count rings in an #LWGEOM.
*/
uint32_t lwgeom_count_rings(const LWGEOM *geom)
{
  int result = 0;

  /* Null? Empty? Zero. */
  if( ! geom || lwgeom_is_empty(geom) )
    return 0;

  switch (geom->type)
  {
  case POINTTYPE:
  case CIRCSTRINGTYPE:
  case COMPOUNDTYPE:
  case MULTICURVETYPE:
  case MULTIPOINTTYPE:
  case MULTILINETYPE:
  case LINETYPE:
    result = 0;
    break;
  case TRIANGLETYPE:
    result = 1;
    break;
  case POLYGONTYPE:
    result = ((LWPOLY *)geom)->nrings;
    break;
  case CURVEPOLYTYPE:
    result = ((LWCURVEPOLY *)geom)->nrings;
    break;
  case MULTISURFACETYPE:
  case MULTIPOLYGONTYPE:
  case POLYHEDRALSURFACETYPE:
  case TINTYPE:
  case COLLECTIONTYPE:
  {
    LWCOLLECTION *col = (LWCOLLECTION*)geom;
    uint32_t i = 0;
    for( i = 0; i < col->ngeoms; i++ )
      result += lwgeom_count_rings(col->geoms[i]);
    break;
  }
  default:
    lwerror("lwgeom_count_rings: unsupported input geometry type: %s", lwtype_name(geom->type));
    break;
  }
  LWDEBUGF(3, "counted %d rings", result);
  return result;
}

int lwgeom_has_srid(const LWGEOM *geom)
{
  if ( geom->srid != SRID_UNKNOWN )
    return LW_TRUE;

  return LW_FALSE;
}


static int lwcollection_dimensionality(const LWCOLLECTION *col)
{
  uint32_t i;
  int dimensionality = 0;
  for ( i = 0; i < col->ngeoms; i++ )
  {
    int d = lwgeom_dimensionality(col->geoms[i]);
    if ( d > dimensionality )
      dimensionality = d;
  }
  return dimensionality;
}

/**
 * Determine the topological dimensionality of a geometry.
 *
 * Returns the maximal topological dimension for the input geometry:
 * - 0 for points and multipoints,
 * - 1 for linear/curve types (lines, circular strings, compound/multi-curves, NURBS),
 * - 2 for polygonal/surface types,
 * Closure-aware dimensionality used by internal validity/repair checks.
 *
 * This is intentionally not the same contract as lwgeom_dimension(), which
 * backs ST_Dimension.  Closed polyhedral surfaces and TINs are reported as 3
 * here so callers that reason about solid-like inputs keep their legacy
 * behavior.
 *
 * For COLLECTIONTYPE the result is the maximum dimensionality of its members.
 *
 * @param geom Geometry to inspect (must be non-NULL).
 * @return The closure-aware dimensionality (0..3). On unsupported types the function logs an error and returns 0.
 */
extern int lwgeom_dimensionality(const LWGEOM *geom)
{
  int dim;

  LWDEBUGF(3, "lwgeom_dimensionality got type %s",
           lwtype_name(geom->type));

  switch (geom->type)
  {
  case POINTTYPE:
  case MULTIPOINTTYPE:
    return 0;
    break;
  case LINETYPE:
  case CIRCSTRINGTYPE:
  case MULTILINETYPE:
  case COMPOUNDTYPE:
  case MULTICURVETYPE:
  case NURBSCURVETYPE:
    return 1;
    break;
  case POLYGONTYPE:
  case TRIANGLETYPE:
  case CURVEPOLYTYPE:
  case MULTIPOLYGONTYPE:
  case MULTISURFACETYPE:
    return 2;
    break;

  case POLYHEDRALSURFACETYPE:
  case TINTYPE:
    dim = lwgeom_is_closed(geom)?3:2;
    return dim;
    break;

  case COLLECTIONTYPE:
    return lwcollection_dimensionality((const LWCOLLECTION *)geom);
    break;
  default:
    lwerror("lwgeom_dimensionality: unsupported input geometry type: %s",
            lwtype_name(geom->type));
    break;
  }
  return 0;
}

extern LWGEOM* lwgeom_remove_repeated_points(const LWGEOM *in, double tolerance)
{
  LWGEOM *out = lwgeom_clone_deep(in);
  lwgeom_remove_repeated_points_in_place(out, tolerance);
  return out;
}

void lwgeom_swap_ordinates(LWGEOM *in, LWORD o1, LWORD o2)
{
  LWCOLLECTION *col;
  LWPOLY *poly;
  uint32_t i;

  if ( (!in) || lwgeom_is_empty(in) ) return;

  /* TODO: check for lwgeom NOT having the specified dimension ? */

  LWDEBUGF(4, "lwgeom_flip_coordinates, got type: %s",
           lwtype_name(in->type));

  switch (in->type)
  {
  case POINTTYPE:
    ptarray_swap_ordinates(lwgeom_as_lwpoint(in)->point, o1, o2);
    break;

  case LINETYPE:
    ptarray_swap_ordinates(lwgeom_as_lwline(in)->points, o1, o2);
    break;

  case CIRCSTRINGTYPE:
    ptarray_swap_ordinates(lwgeom_as_lwcircstring(in)->points, o1, o2);
    break;

  case NURBSCURVETYPE:
  {
    LWNURBSCURVE *n = (LWNURBSCURVE *)in;
    if (n->points)
      ptarray_swap_ordinates(n->points, o1, o2);
    break;
  }

  case POLYGONTYPE:
    poly = (LWPOLY *) in;
    for (i=0; i<poly->nrings; i++)
    {
      ptarray_swap_ordinates(poly->rings[i], o1, o2);
    }
    break;

  case TRIANGLETYPE:
    ptarray_swap_ordinates(lwgeom_as_lwtriangle(in)->points, o1, o2);
    break;

  case MULTIPOINTTYPE:
  case MULTILINETYPE:
  case MULTIPOLYGONTYPE:
  case COLLECTIONTYPE:
  case COMPOUNDTYPE:
  case CURVEPOLYTYPE:
  case MULTISURFACETYPE:
  case MULTICURVETYPE:
  case POLYHEDRALSURFACETYPE:
  case TINTYPE:
    col = (LWCOLLECTION *) in;
    for (i=0; i<col->ngeoms; i++)
    {
      lwgeom_swap_ordinates(col->geoms[i], o1, o2);
    }
    break;

  default:
    lwerror("lwgeom_swap_ordinates: unsupported geometry type: %s",
            lwtype_name(in->type));
    return;
  }

  /* only refresh bbox if X or Y changed */
  if ( in->bbox && (o1 < 2 || o2 < 2) )
  {
    lwgeom_refresh_bbox(in);
  }
}

void lwgeom_set_srid(LWGEOM *geom, int32_t srid)
{
  uint32_t i;

  LWDEBUGF(4,"entered with srid=%d",srid);

  geom->srid = srid;

  if ( lwgeom_is_collection(geom) )
  {
    /* All the children are set to the same SRID value */
    LWCOLLECTION *col = lwgeom_as_lwcollection(geom);
    for ( i = 0; i < col->ngeoms; i++ )
    {
      lwgeom_set_srid(col->geoms[i], srid);
    }
  }
}


/**************************************************************/

static int
cmp_point_x(const void *pa, const void *pb)
{
  LWPOINT *p1 = *((LWPOINT **)pa);
  LWPOINT *p2 = *((LWPOINT **)pb);

  const POINT2D *pt1 = getPoint2d_cp(p1->point, 0);
  const POINT2D *pt2 = getPoint2d_cp(p2->point, 0);

  if (pt1 && pt2)
    return (pt1->x > pt2->x) ? 1 : ((pt1->x < pt2->x) ? -1 : 0);

  if (pt1) return -1;
  if (pt2) return 1;
  return 0;
}

static int
cmp_point_y(const void *pa, const void *pb)
{
  LWPOINT *p1 = *((LWPOINT **)pa);
  LWPOINT *p2 = *((LWPOINT **)pb);

  const POINT2D *pt1 = getPoint2d_cp(p1->point, 0);
  const POINT2D *pt2 = getPoint2d_cp(p2->point, 0);

  if (pt1 && pt2)
    return (pt1->y > pt2->y) ? 1 : ((pt1->y < pt2->y) ? -1 : 0);

  if (pt1) return -1;
  if (pt2) return 1;
  return 0;
}

int
lwgeom_remove_repeated_points_in_place(LWGEOM *geom, double tolerance)
{
  int geometry_modified = LW_FALSE;
  switch (geom->type)
  {
  /* No-op! Cannot remove points */
  case POINTTYPE:
  case TRIANGLETYPE:
    return geometry_modified;
  case LINETYPE: {
    LWLINE *g = (LWLINE *)(geom);
    POINTARRAY *pa = g->points;
    uint32_t npoints = pa->npoints;
    ptarray_remove_repeated_points_in_place(pa, tolerance, 2);
    geometry_modified = npoints != pa->npoints;
    /* Invalid input, discard the collapsed line */
    if (pa->npoints < 2)
    {
      pa->npoints = 0;
      geometry_modified = LW_TRUE;
    }
    break;
  }
  case POLYGONTYPE: {
    uint32_t j = 0;
    LWPOLY *g = (LWPOLY *)(geom);
    for (uint32_t i = 0; i < g->nrings; i++)
    {
      POINTARRAY *pa = g->rings[i];
      uint32_t npoints = pa->npoints;
      ptarray_remove_repeated_points_in_place(pa, tolerance, 4);
      geometry_modified |= npoints != pa->npoints;
      /* Drop collapsed rings */
      if (pa->npoints < 4)
      {
        geometry_modified = LW_TRUE;
        ptarray_free(pa);
        continue;
      }
      g->rings[j++] = pa;
    }
    /* Update ring count */
    g->nrings = j;
    break;
  }
  case MULTIPOINTTYPE: {
    double tolsq = tolerance * tolerance;
    LWMPOINT *mpt = (LWMPOINT *)geom;

    for (uint8_t dim = 0; dim < 2; dim++)
    {
      /* sort by y, then by x - this way the result is sorted by x */
      qsort(mpt->geoms, mpt->ngeoms, sizeof(LWPOINT *), dim ? cmp_point_x : cmp_point_y);
      for (uint32_t i = 0; i < mpt->ngeoms; i++)
      {
        if (!mpt->geoms[i])
          continue;

        const POINT2D *pti = getPoint2d_cp(mpt->geoms[i]->point, 0);
        if (!pti) continue;

        /* check upcoming points if they're within tolerance of current one */
        for (uint32_t j = i + 1; j < mpt->ngeoms; j++)
        {
          if (!mpt->geoms[j])
            continue;

          const POINT2D *ptj = getPoint2d_cp(mpt->geoms[j]->point, 0);
          if (!ptj) continue;

          /* check that the point is in the strip of tolerance around the point */
          if ((dim ? ptj->x - pti->x : ptj->y - pti->y) > tolerance)
            break;

          /* remove any upcoming point that is within tolerance circle */
          if (distance2d_sqr_pt_pt(pti, ptj) <= tolsq)
          {
            lwpoint_free(mpt->geoms[j]);
            mpt->geoms[j] = NULL;
            geometry_modified = LW_TRUE;
          }
        }
      }

      /* null out empties */
      for (uint32_t j = 0; j < mpt->ngeoms; j++)
      {
        if (mpt->geoms[j] && lwpoint_is_empty(mpt->geoms[j]))
        {
          lwpoint_free(mpt->geoms[j]);
          mpt->geoms[j] = NULL;
          geometry_modified = LW_TRUE;
        }
      }

      /* compactify array of points */
      uint32_t i = 0;
      for (uint32_t j = 0; j < mpt->ngeoms; j++)
        if (mpt->geoms[j])
          mpt->geoms[i++] = mpt->geoms[j];
      mpt->ngeoms = i;
    }

    break;
  }

    case CIRCSTRINGTYPE:
    case NURBSCURVETYPE:
      /* Dunno how to handle these, will return untouched */
      return geometry_modified;

  /* Can process most multi* types as generic collection */
  case MULTILINETYPE:
  case MULTIPOLYGONTYPE:
  case TINTYPE:
  case COLLECTIONTYPE:
  /* Curve types we mostly ignore, but allow the linear */
  /* portions to be processed by recursing into them */
  case MULTICURVETYPE:
  case CURVEPOLYTYPE:
  case MULTISURFACETYPE:
  case COMPOUNDTYPE: {
    uint32_t i, j = 0;
    LWCOLLECTION *col = (LWCOLLECTION *)(geom);
    for (i = 0; i < col->ngeoms; i++)
    {
      LWGEOM *g = col->geoms[i];
      if (!g)
        continue;
      geometry_modified |= lwgeom_remove_repeated_points_in_place(g, tolerance);
      /* Drop zero'ed out geometries */
      if (lwgeom_is_empty(g))
      {
        lwgeom_free(g);
        continue;
      }
      col->geoms[j++] = g;
    }
    /* Update geometry count */
    col->ngeoms = j;
    break;
  }
  default: {
    lwerror("%s: unsupported geometry type: %s", __func__, lwtype_name(geom->type));
    break;
  }
  }

  if (geometry_modified)
    lwgeom_drop_bbox(geom);
  return geometry_modified;
}


/**************************************************************/

int
lwgeom_simplify_in_place(LWGEOM *geom, double epsilon, int preserve_collapsed)
{
  int modified = LW_FALSE;

  if (!geom)
  {
    lwerror("NULL geometry encountered");
    return LW_FALSE;
  }

  switch (geom->type)
  {
    /* No-op! Cannot simplify points or triangles */
    case POINTTYPE:
      return modified;
    case TRIANGLETYPE:
    {
      if (preserve_collapsed)
        return modified;
      LWTRIANGLE *t = lwgeom_as_lwtriangle(geom);
      POINTARRAY *pa = t->points;
      ptarray_simplify_in_place(pa, epsilon, 0);
      if (pa->npoints < 3)
      {
        pa->npoints = 0;
        modified = LW_TRUE;
      }
      break;
    }
    case LINETYPE:
    {
      LWLINE *g = (LWLINE*)(geom);
      POINTARRAY *pa = g->points;
      uint32_t in_npoints = pa->npoints;
      ptarray_simplify_in_place(pa, epsilon, 2);
      modified = in_npoints != pa->npoints;
      /* Invalid output */
      if (pa->npoints == 1 && pa->maxpoints > 1)
      {
        /* Use first point as last point */
        if (preserve_collapsed)
        {
          pa->npoints = 2;
          ptarray_copy_point(pa, 0, 1);
        }
        /* Finish the collapse process */
        else
        {
          pa->npoints = 0;
        }
      }
      /* Duped output, force collapse */
      if (pa->npoints == 2 && !preserve_collapsed)
      {
        if (p2d_same(getPoint2d_cp(pa, 0), getPoint2d_cp(pa, 1)))
          pa->npoints = 0;
      }
      break;
    }
    case POLYGONTYPE:
    {
      uint32_t i, j = 0;
      LWPOLY *g = (LWPOLY*)(geom);
      for (i = 0; i < g->nrings; i++)
      {
        POINTARRAY *pa = g->rings[i];
        /* Only stop collapse on first ring */
        int minpoints = (preserve_collapsed && i == 0) ? 4 : 0;
        /* Skip zero'ed out rings */
        if(!pa)
          continue;
        uint32_t in_npoints = pa->npoints;
        ptarray_simplify_in_place(pa, epsilon, minpoints);
        modified |= in_npoints != pa->npoints;
        /* Drop collapsed rings */
        if(pa->npoints < 4)
        {
          if (i == 0)
          {
            /* If the outer ring is dropped, all can be dropped */
            for (i = 0; i < g->nrings; i++)
            {
              pa = g->rings[i];
              ptarray_free(pa);
            }
            break;
          }
          else
          {
            /* Drop this inner ring only */
            ptarray_free(pa);
            continue;
          }
        }
        g->rings[j++] = pa;
      }
      /* Update ring count */
      g->nrings = j;
      break;
    }
    /* Can process all multi* types as generic collection */
    case MULTIPOINTTYPE:
    case MULTILINETYPE:
    case MULTIPOLYGONTYPE:
    case TINTYPE:
    case COLLECTIONTYPE:
    {
      uint32_t i, j = 0;
      LWCOLLECTION *col = (LWCOLLECTION*)geom;
      for (i = 0; i < col->ngeoms; i++)
      {
        LWGEOM *g = col->geoms[i];
        if (!g) continue;
        modified |= lwgeom_simplify_in_place(g, epsilon, preserve_collapsed);
        /* Drop zero'ed out geometries */
        if(lwgeom_is_empty(g))
        {
          lwgeom_free(g);
          continue;
        }
        col->geoms[j++] = g;
      }
      /* Update geometry count */
      col->ngeoms = j;
      break;
    }
    default:
    {
      lwerror("%s: unsupported geometry type: %s", __func__, lwtype_name(geom->type));
      break;
    }
  }

  if (modified)
  {
    lwgeom_drop_bbox(geom);
  }
  return modified;
}

LWGEOM* lwgeom_simplify(const LWGEOM *igeom, double dist, int preserve_collapsed)
{
  LWGEOM *lwgeom_out = lwgeom_clone_deep(igeom);
  lwgeom_simplify_in_place(lwgeom_out, dist, preserve_collapsed);
  if (lwgeom_is_empty(lwgeom_out))
  {
    lwgeom_free(lwgeom_out);
    return NULL;
  }
  return lwgeom_out;
}

/**************************************************************/


double lwgeom_area(const LWGEOM *geom)
{
  int type = geom->type;

  if ( type == POLYGONTYPE )
    return lwpoly_area((LWPOLY*)geom);
  else if ( type == CURVEPOLYTYPE )
    return lwcurvepoly_area((LWCURVEPOLY*)geom);
  else if (type ==  TRIANGLETYPE )
    return lwtriangle_area((LWTRIANGLE*)geom);
  else if ( lwgeom_is_collection(geom) )
  {
    double area = 0.0;
    uint32_t i;
    LWCOLLECTION *col = (LWCOLLECTION*)geom;
    for ( i = 0; i < col->ngeoms; i++ )
      area += lwgeom_area(col->geoms[i]);
    return area;
  }
  else
    return 0.0;
}

double lwgeom_perimeter(const LWGEOM *geom)
{
  int type = geom->type;
  if ( type == POLYGONTYPE )
    return lwpoly_perimeter((LWPOLY*)geom);
  else if ( type == CURVEPOLYTYPE )
    return lwcurvepoly_perimeter((LWCURVEPOLY*)geom);
  else if ( type == TRIANGLETYPE )
    return lwtriangle_perimeter((LWTRIANGLE*)geom);
  else if ( lwgeom_is_collection(geom) )
  {
    double perimeter = 0.0;
    uint32_t i;
    LWCOLLECTION *col = (LWCOLLECTION*)geom;
    for ( i = 0; i < col->ngeoms; i++ )
      perimeter += lwgeom_perimeter(col->geoms[i]);
    return perimeter;
  }
  else
    return 0.0;
}

double lwgeom_perimeter_2d(const LWGEOM *geom)
{
  int type = geom->type;
  if ( type == POLYGONTYPE )
    return lwpoly_perimeter_2d((LWPOLY*)geom);
  else if ( type == CURVEPOLYTYPE )
    return lwcurvepoly_perimeter_2d((LWCURVEPOLY*)geom);
  else if ( type == TRIANGLETYPE )
    return lwtriangle_perimeter_2d((LWTRIANGLE*)geom);
  else if ( lwgeom_is_collection(geom) )
  {
    double perimeter = 0.0;
    uint32_t i;
    LWCOLLECTION *col = (LWCOLLECTION*)geom;
    for ( i = 0; i < col->ngeoms; i++ )
      perimeter += lwgeom_perimeter_2d(col->geoms[i]);
    return perimeter;
  }
  else
    return 0.0;
}

double lwgeom_length(const LWGEOM *geom)
{
  int type = geom->type;
  if ( type == LINETYPE )
    return lwline_length((LWLINE*)geom);
  else if ( type == CIRCSTRINGTYPE )
    return lwcircstring_length((LWCIRCSTRING*)geom);
  else if ( type == NURBSCURVETYPE )
  {
    double length;
    LWLINE *line = lwnurbscurve_to_linestring((LWNURBSCURVE*)geom, 32);
    length = lwline_length(line);
    lwline_free(line);
    return length;
  }
  else if ( type == COMPOUNDTYPE )
    return lwcompound_length((LWCOMPOUND*)geom);
  else if ( lwgeom_is_collection(geom) )
  {
    double length = 0.0;
    uint32_t i;
    LWCOLLECTION *col = (LWCOLLECTION*)geom;
    for ( i = 0; i < col->ngeoms; i++ )
      length += lwgeom_length(col->geoms[i]);
    return length;
  }
  else
    return 0.0;
}

double lwgeom_length_2d(const LWGEOM *geom)
{
  int type = geom->type;
  if ( type == LINETYPE )
    return lwline_length_2d((LWLINE*)geom);
  else if ( type == CIRCSTRINGTYPE )
    return lwcircstring_length_2d((LWCIRCSTRING*)geom);
  else if ( type == NURBSCURVETYPE )
  {
    double length;
    LWLINE *line = lwnurbscurve_to_linestring((LWNURBSCURVE*)geom, 32);
    length = lwline_length_2d(line);
    lwline_free(line);
    return length;
  }
  else if ( type == COMPOUNDTYPE )
    return lwcompound_length_2d((LWCOMPOUND*)geom);
  else if ( type != CURVEPOLYTYPE && lwgeom_is_collection(geom) )
  {
    double length = 0.0;
    uint32_t i;
    LWCOLLECTION *col = (LWCOLLECTION*)geom;
    for ( i = 0; i < col->ngeoms; i++ )
      length += lwgeom_length_2d(col->geoms[i]);
    return length;
  }
  else
    return 0.0;
}

void
lwgeom_affine(LWGEOM *geom, const AFFINE *affine)
{
  int type = geom->type;
  uint32_t i;

  switch(type)
  {
    /* Take advantage of fact that pt/ln/circ/tri have same memory structure */
    case POINTTYPE:
    case LINETYPE:
    case CIRCSTRINGTYPE:
    case TRIANGLETYPE:
    {
      LWLINE *l = (LWLINE*)geom;
      ptarray_affine(l->points, affine);
      break;
    }
    case NURBSCURVETYPE:
    {
      LWNURBSCURVE *n = (LWNURBSCURVE*)geom;
      if (n->points != NULL)
        ptarray_affine(n->points, affine);
      break;
    }
    case POLYGONTYPE:
    {
      LWPOLY *p = (LWPOLY*)geom;
      for( i = 0; i < p->nrings; i++ )
        ptarray_affine(p->rings[i], affine);
      break;
    }
    case CURVEPOLYTYPE:
    {
      LWCURVEPOLY *c = (LWCURVEPOLY*)geom;
      for( i = 0; i < c->nrings; i++ )
        lwgeom_affine(c->rings[i], affine);
      break;
    }
    default:
    {
      if( lwgeom_is_collection(geom) )
      {
        LWCOLLECTION *c = (LWCOLLECTION*)geom;
        for( i = 0; i < c->ngeoms; i++ )
        {
          lwgeom_affine(c->geoms[i], affine);
        }
      }
      else
      {
        lwerror("lwgeom_affine: unable to handle type '%s'", lwtype_name(type));
      }
    }
  }

  /* Recompute bbox if needed */
  if (geom->bbox)
    lwgeom_refresh_bbox(geom);
}

/**
 * Scale a geometry in-place by per-ordinate factors.
 *
 * Applies the provided scale factors to all coordinates of the given geometry.
 * Scaling is applied component-wise: X by factor->x, Y by factor->y, Z by factor->z
 * and M by factor->m when those ordinates are present. The operation mutates
 * the input geometry and recurses into collections and curve/polygon components.
 * If the geometry has an attached bounding box it is refreshed after scaling.
 *
 * @param geom Geometry to be scaled (modified in place). Behavior is undefined if NULL.
 * @param factor Per-ordinate scale factors (X, Y, Z, M).
 */
void
lwgeom_scale(LWGEOM *geom, const POINT4D *factor)
{
  int type = geom->type;
  uint32_t i;

  switch(type)
  {
    /* Take advantage of fact that pt/ln/circ/tri have same memory structure */
    case POINTTYPE:
    case LINETYPE:
    case CIRCSTRINGTYPE:
    case TRIANGLETYPE:
    {
      LWLINE *l = (LWLINE*)geom;
      ptarray_scale(l->points, factor);
      break;
    }
    case NURBSCURVETYPE:
    {
      LWNURBSCURVE *n = (LWNURBSCURVE*)geom;
      if (n->points != NULL)
        ptarray_scale(n->points, factor);
      break;
    }
    case POLYGONTYPE:
    {
      LWPOLY *p = (LWPOLY*)geom;
      for( i = 0; i < p->nrings; i++ )
        ptarray_scale(p->rings[i], factor);
      break;
    }
    case CURVEPOLYTYPE:
    {
      LWCURVEPOLY *c = (LWCURVEPOLY*)geom;
      for( i = 0; i < c->nrings; i++ )
        lwgeom_scale(c->rings[i], factor);
      break;
    }
    default:
    {
      if( lwgeom_is_collection(geom) )
      {
        LWCOLLECTION *c = (LWCOLLECTION*)geom;
        for( i = 0; i < c->ngeoms; i++ )
        {
          lwgeom_scale(c->geoms[i], factor);
        }
      }
      else
      {
        lwerror("lwgeom_scale: unable to handle type '%s'", lwtype_name(type));
      }
    }
  }

  /* Recompute bbox if needed */
  if (geom->bbox)
    lwgeom_refresh_bbox(geom);
}

/**
 * Construct an empty LWGEOM of the requested geometry type.
 *
 * Creates and returns a newly allocated geometry instance with no coordinates
 * (an "empty" geometry), using the provided SRID and dimensionality flags.
 * For homogeneous collection types (MULTI*, COLLECTION, COMPOUND) the `type`
 * argument selects the specific collection subtype to construct.
 *
 * @param type Geometry type code (e.g., POINTTYPE, LINETYPE, POLYGONTYPE,
 *             CURVEPOLYTYPE, CIRCSTRINGTYPE, TRIANGLETYPE, COMPOUNDTYPE,
 *             MULTIPOINTTYPE, MULTILINETYPE, MULTIPOLYGONTYPE, COLLECTIONTYPE,
 *             NURBSCURVETYPE). Unsupported types cause an error to be logged
 *             and the function to return NULL.
 * @param srid Spatial reference identifier to assign to the constructed geometry
 *             (use SRID_UNKNOWN if not set).
 * @param hasz Non-zero to include a Z ordinate in the geometry, zero to omit Z.
 * @param hasm Non-zero to include an M ordinate in the geometry, zero to omit M.
 *
 * @return Pointer to the newly constructed LWGEOM, or NULL if `type` is not
 *         supported.
 */
LWGEOM *
lwgeom_construct_empty(uint8_t type, int32_t srid, char hasz, char hasm)
{
  switch(type)
  {
    case POINTTYPE:
      return lwpoint_as_lwgeom(lwpoint_construct_empty(srid, hasz, hasm));
    case LINETYPE:
      return lwline_as_lwgeom(lwline_construct_empty(srid, hasz, hasm));
    case POLYGONTYPE:
      return lwpoly_as_lwgeom(lwpoly_construct_empty(srid, hasz, hasm));
    case CURVEPOLYTYPE:
      return lwcurvepoly_as_lwgeom(lwcurvepoly_construct_empty(srid, hasz, hasm));
    case CIRCSTRINGTYPE:
      return lwcircstring_as_lwgeom(lwcircstring_construct_empty(srid, hasz, hasm));
    case TRIANGLETYPE:
      return lwtriangle_as_lwgeom(lwtriangle_construct_empty(srid, hasz, hasm));
    case COMPOUNDTYPE:
    case MULTIPOINTTYPE:
    case MULTILINETYPE:
    case MULTIPOLYGONTYPE:
    case COLLECTIONTYPE:
      return lwcollection_as_lwgeom(lwcollection_construct_empty(type, srid, hasz, hasm));
    case NURBSCURVETYPE:
      return lwnurbscurve_as_lwgeom(lwnurbscurve_construct_empty(srid, hasz, hasm));
    default:
      lwerror("lwgeom_construct_empty: unsupported geometry type: %s",
              lwtype_name(type));
      return NULL;
  }
}

int
lwgeom_startpoint(const LWGEOM *lwgeom, POINT4D *pt)
{
  if ( ! lwgeom || lwgeom_is_empty(lwgeom) )
    return LW_FAILURE;

  switch( lwgeom->type )
  {
    case POINTTYPE:
      return ptarray_startpoint(((LWPOINT*)lwgeom)->point, pt);
    case TRIANGLETYPE:
    case CIRCSTRINGTYPE:
    case LINETYPE:
      return ptarray_startpoint(((LWLINE*)lwgeom)->points, pt);
    case NURBSCURVETYPE:
    {
      LWPOINT *start = lwnurbscurve_evaluate((const LWNURBSCURVE*)lwgeom, 0.0);
      int rv = LW_FAILURE;
      if ( start && start->point && start->point->npoints > 0 )
        rv = getPoint4d_p(start->point, 0, pt) ? LW_SUCCESS : LW_FAILURE;
      lwpoint_free(start);
      return rv;
    }
    case POLYGONTYPE:
      return lwpoly_startpoint((LWPOLY*)lwgeom, pt);
    case TINTYPE:
    case CURVEPOLYTYPE:
    case COMPOUNDTYPE:
    case MULTIPOINTTYPE:
    case MULTILINETYPE:
    case MULTIPOLYGONTYPE:
    case COLLECTIONTYPE:
    case POLYHEDRALSURFACETYPE:
      return lwcollection_startpoint((LWCOLLECTION*)lwgeom, pt);
    default:
      lwerror("lwgeom_startpoint: unsupported geometry type: %s", lwtype_name(lwgeom->type));
      return LW_FAILURE;
  }
}

static inline double
snap_to_int(double val)
{
  const double tolerance = 1e-6;
    double rintval = rint(val);
    if (fabs(val - rintval) < tolerance)
    {
        return rintval;
    }
    return val;
}

/*
 * See https://github.com/libgeos/geos/pull/956
 * We use scale for rounding when gridsize is < 1 and
 * gridsize for rounding when scale < 1.
 */
static inline void
condition_gridspec_scale(gridspec *grid)
{
  if(grid->xsize > 0)
  {
    if(grid->xsize < 1)
      grid->xscale = snap_to_int(1/grid->xsize);
    else
      grid->xsize = snap_to_int(grid->xsize);
  }
  if(grid->ysize > 0)
  {
    if(grid->ysize < 1)
      grid->yscale = snap_to_int(1/grid->ysize);
    else
      grid->ysize = snap_to_int(grid->ysize);
  }
}


void
lwgeom_grid_in_place(LWGEOM *geom, gridspec *grid)
{
  if (!geom) return;
  if (lwgeom_is_empty(geom)) return;

  condition_gridspec_scale(grid);

  switch ( geom->type )
  {
    case POINTTYPE:
    {
      LWPOINT *pt = (LWPOINT*)(geom);
      ptarray_grid_in_place(pt->point, grid);
      return;
    }
    case CIRCSTRINGTYPE:
    case TRIANGLETYPE:
    case LINETYPE:
    {
      LWLINE *ln = (LWLINE*)(geom);
      ptarray_grid_in_place(ln->points, grid);
      /* For invalid line, return an EMPTY */
      if (ln->points->npoints < 2)
        ln->points->npoints = 0;
      return;
    }
    case NURBSCURVETYPE:
    {
      LWNURBSCURVE *n = (LWNURBSCURVE *)(geom);
      uint32_t old_npoints;
      if (!n->points)
        return;

      old_npoints = n->points->npoints;
      ptarray_grid_in_place(n->points, grid);
      if (n->points->npoints < n->degree + 1)
        n->points->npoints = 0;

      if (n->points->npoints != old_npoints)
      {
        if (n->weights)
          lwfree(n->weights);
        if (n->knots)
          lwfree(n->knots);
        n->weights = NULL;
        n->knots = NULL;
        n->nweights = 0;
        n->nknots = 0;
      }
      return;
    }
    case POLYGONTYPE:
    {
      LWPOLY *ply = (LWPOLY*)(geom);
      if (!ply->rings) return;

      /* Check first the external ring */
      uint32_t i = 0;
      POINTARRAY *pa = ply->rings[0];
      ptarray_grid_in_place(pa, grid);
      if (pa->npoints < 4)
      {
        /* External ring collapsed: free everything */
        for (i = 0; i < ply->nrings; i++)
        {
          ptarray_free(ply->rings[i]);
        }
        ply->nrings = 0;
        return;
      }

      /* Check the other rings */
      uint32_t j = 1;
      for (i = 1; i < ply->nrings; i++)
      {
        POINTARRAY *pa = ply->rings[i];
        ptarray_grid_in_place(pa, grid);

        /* Skip bad rings */
        if (pa->npoints >= 4)
        {
          ply->rings[j++] = pa;
        }
        else
        {
          ptarray_free(pa);
        }
      }
      /* Adjust ring count appropriately */
      ply->nrings = j;
      return;
    }
    case MULTIPOINTTYPE:
    case MULTILINETYPE:
    case MULTIPOLYGONTYPE:
    case TINTYPE:
    case COLLECTIONTYPE:
    case COMPOUNDTYPE:
    {
      LWCOLLECTION *col = (LWCOLLECTION*)(geom);
      uint32_t i, j = 0;
      if (!col->geoms) return;
      for (i = 0; i < col->ngeoms; i++)
      {
        LWGEOM *g = col->geoms[i];
        lwgeom_grid_in_place(g, grid);
        /* Empty geoms need to be freed */
        /* before we move on */
        if (lwgeom_is_empty(g))
        {
          lwgeom_free(g);
          continue;
        }
        col->geoms[j++] = g;
      }
      col->ngeoms = j;
      return;
    }
    default:
    {
      lwerror("%s: Unsupported geometry type: %s", __func__,
              lwtype_name(geom->type));
      return;
    }
  }
}


LWGEOM *
lwgeom_grid(const LWGEOM *lwgeom, gridspec *grid)
{
  LWGEOM *lwgeom_out = lwgeom_clone_deep(lwgeom);
  lwgeom_grid_in_place(lwgeom_out, grid);
  return lwgeom_out;
}


/* Prototype for recursion */
static void lwgeom_subdivide_recursive(const LWGEOM *geom,
               uint8_t dimension,
               uint32_t maxvertices,
               uint32_t depth,
               LWCOLLECTION *col,
               double gridSize);

static void
lwgeom_subdivide_recursive(const LWGEOM *geom,
         uint8_t dimension,
         uint32_t maxvertices,
         uint32_t depth,
         LWCOLLECTION *col,
         double gridSize)
{
  const uint32_t maxdepth = 50;

  if (!geom)
    return;

  const GBOX *box_in = lwgeom_get_bbox(geom);
  if (!box_in)
    return;

  LW_ON_INTERRUPT(return);

  GBOX clip;
  gbox_duplicate(box_in, &clip);
  double width = clip.xmax - clip.xmin;
  double height = clip.ymax - clip.ymin;

  if ( geom->type == POLYHEDRALSURFACETYPE || geom->type == TINTYPE )
    lwerror("%s: unsupported geometry type '%s'", __func__, lwtype_name(geom->type));

  if ( width == 0.0 && height == 0.0 )
  {
    if ( geom->type == POINTTYPE && dimension == 0)
      lwcollection_add_lwgeom(col, lwgeom_clone_deep(geom));
    return;
  }

  if (width == 0.0)
  {
    clip.xmax += FP_TOLERANCE;
    clip.xmin -= FP_TOLERANCE;
    width = 2 * FP_TOLERANCE;
  }
  if (height == 0.0)
  {
    clip.ymax += FP_TOLERANCE;
    clip.ymin -= FP_TOLERANCE;
    height = 2 * FP_TOLERANCE;
  }

  /* Always just recurse into collections */
  if ( lwgeom_is_collection(geom) && geom->type != MULTIPOINTTYPE )
  {
    LWCOLLECTION *incol = (LWCOLLECTION*)geom;
    /* Don't increment depth yet, since we aren't actually
     * subdividing geometries yet */
    for (uint32_t i = 0; i < incol->ngeoms; i++ )
      lwgeom_subdivide_recursive(incol->geoms[i], dimension, maxvertices, depth, col, gridSize);
    return;
  }

  if (lwgeom_dimension(geom) < dimension)
  {
    /* We've hit a lower dimension object produced by clipping at
     * a shallower recursion level. Ignore it. */
    return;
  }

  /* But don't go too far. 2^50 ~= 10^15, that's enough subdivision */
  /* Just add what's left */
  if ( depth > maxdepth )
  {
    lwcollection_add_lwgeom(col, lwgeom_clone_deep(geom));
    return;
  }

  uint32_t nvertices = lwgeom_count_vertices(geom);

  /* Skip empties entirely */
  if (nvertices == 0)
    return;

  /* If it is under the vertex tolerance, just add it, we're done */
  if (nvertices <= maxvertices)
  {
    lwcollection_add_lwgeom(col, lwgeom_clone_deep(geom));
    return;
  }

  uint8_t split_ordinate = (width > height) ? 0 : 1;
  double center = (split_ordinate == 0) ? (clip.xmin + clip.xmax) / 2 : (clip.ymin + clip.ymax) / 2;
  double pivot = DBL_MAX;
  if (geom->type == POLYGONTYPE)
  {
    uint32_t ring_to_trim = 0;
    double ring_area = 0;
    double pivot_eps = DBL_MAX;
    double pt_eps = DBL_MAX;
    POINTARRAY *pa;
    LWPOLY *lwpoly = (LWPOLY *)geom;

    /* if there are more points in holes than in outer ring */
    if (nvertices >= 2 * lwpoly->rings[0]->npoints)
    {
      /* trim holes starting from biggest */
      for (uint32_t i = 1; i < lwpoly->nrings; i++)
      {
        double current_ring_area = fabs(ptarray_signed_area(lwpoly->rings[i]));
        if (current_ring_area >= ring_area)
        {
          ring_area = current_ring_area;
          ring_to_trim = i;
        }
      }
    }

    pa = lwpoly->rings[ring_to_trim];

    /* find most central point in chosen ring */
    for (uint32_t i = 0; i < pa->npoints; i++)
    {
      double pt;
      if (split_ordinate == 0)
        pt = getPoint2d_cp(pa, i)->x;
      else
        pt = getPoint2d_cp(pa, i)->y;
      pt_eps = fabs(pt - center);
      if (pivot_eps > pt_eps)
      {
        pivot = pt;
        pivot_eps = pt_eps;
      }
    }
  }
  GBOX subbox1, subbox2;
  gbox_duplicate(&clip, &subbox1);
  gbox_duplicate(&clip, &subbox2);

  if (pivot == DBL_MAX)
    pivot = center;

  if (split_ordinate == 0)
  {
    if (FP_NEQUALS(subbox1.xmax, pivot) && FP_NEQUALS(subbox1.xmin, pivot))
      subbox1.xmax = subbox2.xmin = pivot;
    else
      subbox1.xmax = subbox2.xmin = center;
  }
  else
  {
    if (FP_NEQUALS(subbox1.ymax, pivot) && FP_NEQUALS(subbox1.ymin, pivot))
      subbox1.ymax = subbox2.ymin = pivot;
    else
      subbox1.ymax = subbox2.ymin = center;
  }

  ++depth;

  {
    LWGEOM *subbox = (LWGEOM *)lwpoly_construct_envelope(
        geom->srid, subbox1.xmin, subbox1.ymin, subbox1.xmax, subbox1.ymax);
    LWGEOM *clipped = lwgeom_intersection_prec(geom, subbox, gridSize);
    lwgeom_simplify_in_place(clipped, 0.0, LW_TRUE);
    lwgeom_free(subbox);
    if (clipped && !lwgeom_is_empty(clipped))
    {
      lwgeom_subdivide_recursive(clipped, dimension, maxvertices, depth, col, gridSize);
      lwgeom_free(clipped);
    }
  }
  {
    LWGEOM *subbox = (LWGEOM *)lwpoly_construct_envelope(
        geom->srid, subbox2.xmin, subbox2.ymin, subbox2.xmax, subbox2.ymax);
    LWGEOM *clipped = lwgeom_intersection_prec(geom, subbox, gridSize);
    lwgeom_simplify_in_place(clipped, 0.0, LW_TRUE);
    lwgeom_free(subbox);
    if (clipped && !lwgeom_is_empty(clipped))
    {
      lwgeom_subdivide_recursive(clipped, dimension, maxvertices, depth, col, gridSize);
      lwgeom_free(clipped);
    }
  }
}

LWCOLLECTION *
lwgeom_subdivide_prec(const LWGEOM *geom, uint32_t maxvertices, double gridSize)
{
  static uint32_t startdepth = 0;
  static uint32_t minmaxvertices = 5;
  LWCOLLECTION *col;

  col = lwcollection_construct_empty(COLLECTIONTYPE, geom->srid, lwgeom_has_z(geom), lwgeom_has_m(geom));

  if ( lwgeom_is_empty(geom) )
    return col;

  if ( maxvertices < minmaxvertices )
  {
    lwcollection_free(col);
    lwerror("%s: cannot subdivide to fewer than %d vertices per output", __func__, minmaxvertices);
  }

  lwgeom_subdivide_recursive(geom, lwgeom_dimension(geom), maxvertices, startdepth, col, gridSize);
  lwgeom_set_srid((LWGEOM*)col, geom->srid);
  return col;
}

LWCOLLECTION *
lwgeom_subdivide(const LWGEOM *geom, uint32_t maxvertices)
{
  return lwgeom_subdivide_prec(geom, maxvertices, -1);
}



int
lwgeom_is_trajectory(const LWGEOM *geom)
{
  int type = geom->type;

  if( type != LINETYPE )
  {
    lwnotice("Geometry is not a LINESTRING");
    return LW_FALSE;
  }
  return lwline_is_trajectory((LWLINE*)geom);
}

#define STATIC_ASSERT(COND,MSG) typedef char static_assertion_##MSG[(COND)?1:-1]
STATIC_ASSERT(sizeof(double) == sizeof(uint64_t),this_should_be_true);

static double trim_preserve_decimal_digits(double d, int32_t decimal_digits)
{
  uint64_t dint = 0;
  memcpy(&dint, &d, sizeof(double));
  /* Extract the exponent from the IEEE 754 integer representation, which */
  /* corresponds to floor(log2(fabs(d))) */
  const int exponent = (int)((dint >> 52) & 2047) - 1023;
  /* (x * 851 + 255) / 256 == 1 + (int)(x * log2(10)) for x in [0,30] */
  int bits_needed = 1 + exponent + (decimal_digits * 851 + 255) / 256;
  /* for negative values, (x * 851 + 255) / 256 == (int)(x * log2(10)), so */
  /* subtract one */
  if (decimal_digits < 0)
    bits_needed --;

  /* This will also handle NaN and Inf since exponent = 1023, and thus for */
  /* reasonable decimal_digits values bits_needed will be > 52 */
  if (bits_needed >= 52)
  {
    return d;
  }
  if (bits_needed < 1 )
    bits_needed = 1;
  const uint64_t mask = 0xffffffffffffffffULL << (52 - bits_needed);
  dint &= mask;
  memcpy(&d, &dint, sizeof(double));
  return d;
}

void lwgeom_trim_bits_in_place(LWGEOM* geom, int32_t prec_x, int32_t prec_y, int32_t prec_z, int32_t prec_m)
{
  LWPOINTITERATOR* it = lwpointiterator_create_rw(geom);
  POINT4D p;

  while (lwpointiterator_has_next(it))
  {
    lwpointiterator_peek(it, &p);
    p.x = trim_preserve_decimal_digits(p.x, prec_x);
    p.y = trim_preserve_decimal_digits(p.y, prec_y);
    if (lwgeom_has_z(geom))
      p.z = trim_preserve_decimal_digits(p.z, prec_z);
    if (lwgeom_has_m(geom))
      p.m = trim_preserve_decimal_digits(p.m, prec_m);
    lwpointiterator_modify_next(it, &p);
  }

  lwpointiterator_destroy(it);
}

LWGEOM *
lwgeom_boundary(LWGEOM *lwgeom)
{
  int32_t srid = lwgeom_get_srid(lwgeom);
  uint8_t hasz = lwgeom_has_z(lwgeom);
  uint8_t hasm = lwgeom_has_m(lwgeom);

  switch (lwgeom->type)
  {
  case POINTTYPE:
  case MULTIPOINTTYPE: {
    return lwgeom_construct_empty(lwgeom->type, srid, hasz, hasm);
  }
  case LINETYPE:
  case CIRCSTRINGTYPE: {
    if (lwgeom_is_closed(lwgeom) || lwgeom_is_empty(lwgeom))
      return (LWGEOM *)lwmpoint_construct_empty(srid, hasz, hasm);
    else
    {
      LWLINE *lwline = (LWLINE *)lwgeom;
      LWMPOINT *lwmpoint = lwmpoint_construct_empty(srid, hasz, hasm);
      POINT4D pt;
      getPoint4d_p(lwline->points, 0, &pt);
      lwmpoint_add_lwpoint(lwmpoint, lwpoint_make(srid, hasz, hasm, &pt));
      getPoint4d_p(lwline->points, lwline->points->npoints - 1, &pt);
      lwmpoint_add_lwpoint(lwmpoint, lwpoint_make(srid, hasz, hasm, &pt));

      return (LWGEOM *)lwmpoint;
    }
  }
  case NURBSCURVETYPE: {
    if (lwgeom_is_closed(lwgeom) || lwgeom_is_empty(lwgeom))
      return (LWGEOM *)lwmpoint_construct_empty(srid, hasz, hasm);
    else
    {
      LWMPOINT *lwmpoint = lwmpoint_construct_empty(srid, hasz, hasm);
      POINT4D pt;
      LWPOINT *startpt = lwnurbscurve_evaluate((const LWNURBSCURVE *)lwgeom, 0.0);
      LWPOINT *endpt   = lwnurbscurve_evaluate((const LWNURBSCURVE *)lwgeom, 1.0);
      if (startpt)
      {
        lwpoint_getPoint4d_p(startpt, &pt);
        lwmpoint_add_lwpoint(lwmpoint, lwpoint_make(srid, hasz, hasm, &pt));
        lwpoint_free(startpt);
      }
      if (endpt)
      {
        lwpoint_getPoint4d_p(endpt, &pt);
        lwmpoint_add_lwpoint(lwmpoint, lwpoint_make(srid, hasz, hasm, &pt));
        lwpoint_free(endpt);
      }
      return (LWGEOM *)lwmpoint;
    }
  }
  case MULTILINETYPE:
  case MULTICURVETYPE: {
    LWMLINE *lwmline = (LWMLINE *)lwgeom;
    POINT4D *out = lwalloc(sizeof(POINT4D) * lwmline->ngeoms * 2);
    uint32_t n = 0;

    for (uint32_t i = 0; i < lwmline->ngeoms; i++)
    {
      LWMPOINT *points = lwgeom_as_lwmpoint(lwgeom_boundary((LWGEOM *)lwmline->geoms[i]));
      if (!points)
        continue;

      for (uint32_t k = 0; k < points->ngeoms; k++)
      {
        POINT4D pt = getPoint4d(points->geoms[k]->point, 0);

        uint8_t seen = LW_FALSE;
        for (uint32_t j = 0; j < n; j++)
        {
          if (memcmp(&(out[j]), &pt, sizeof(POINT4D)) == 0)
          {
            seen = LW_TRUE;
            out[j] = out[--n];
            break;
          }
        }
        if (!seen)
          out[n++] = pt;
      }

      lwgeom_free((LWGEOM *)points);
    }

    LWMPOINT *lwmpoint = lwmpoint_construct_empty(srid, hasz, hasm);

    for (uint32_t i = 0; i < n; i++)
      lwmpoint_add_lwpoint(lwmpoint, lwpoint_make(srid, hasz, hasm, &(out[i])));

    lwfree(out);

    return (LWGEOM *)lwmpoint;
  }
  case TRIANGLETYPE: {
    LWTRIANGLE *lwtriangle = (LWTRIANGLE *)lwgeom;
    POINTARRAY *points = ptarray_clone_deep(lwtriangle->points);
    return (LWGEOM *)lwline_construct(srid, 0, points);
  }
  case POLYGONTYPE: {
    LWPOLY *lwpoly = (LWPOLY *)lwgeom;

    LWMLINE *lwmline = lwmline_construct_empty(srid, hasz, hasm);
    for (uint32_t i = 0; i < lwpoly->nrings; i++)
    {
      POINTARRAY *ring = ptarray_clone_deep(lwpoly->rings[i]);
      lwmline_add_lwline(lwmline, lwline_construct(srid, 0, ring));
    }

    /* Homogenize the multilinestring to hopefully get a single LINESTRING */
    LWGEOM *lwout = lwgeom_homogenize((LWGEOM *)lwmline);
    lwgeom_free((LWGEOM *)lwmline);
    return lwout;
  }
  case CURVEPOLYTYPE: {
    LWCURVEPOLY *lwcurvepoly = (LWCURVEPOLY *)lwgeom;
    LWCOLLECTION *lwcol = lwcollection_construct_empty(MULTICURVETYPE, srid, hasz, hasm);

    for (uint32_t i = 0; i < lwcurvepoly->nrings; i++)
      lwcol = lwcollection_add_lwgeom(lwcol, lwgeom_clone_deep(lwcurvepoly->rings[i]));

    return (LWGEOM *)lwcol;
  }
  case MULTIPOLYGONTYPE:
  case COLLECTIONTYPE:
  case TINTYPE: {
    LWCOLLECTION *lwcol = (LWCOLLECTION *)lwgeom;
    LWCOLLECTION *lwcol_boundary = lwcollection_construct_empty(COLLECTIONTYPE, srid, hasz, hasm);

    for (uint32_t i = 0; i < lwcol->ngeoms; i++)
      lwcollection_add_lwgeom(lwcol_boundary, lwgeom_boundary(lwcol->geoms[i]));

    LWGEOM *lwout = lwgeom_homogenize((LWGEOM *)lwcol_boundary);
    lwgeom_free((LWGEOM *)lwcol_boundary);

    return lwout;
  }
  default:
    lwerror("%s: unsupported geometry type: %s", __func__, lwtype_name(lwgeom->type));
    return NULL;
  }
}

int
lwgeom_isfinite(const LWGEOM *lwgeom)
{
  LWPOINTITERATOR* it = lwpointiterator_create(lwgeom);
  int hasz = lwgeom_has_z(lwgeom);
  int hasm = lwgeom_has_m(lwgeom);

  while (lwpointiterator_has_next(it))
  {
    POINT4D p;
    lwpointiterator_next(it, &p);
    int finite = isfinite(p.x) &&
      isfinite(p.y) &&
      (hasz ? isfinite(p.z) : 1) &&
      (hasm ? isfinite(p.m) : 1);

    if (!finite)
    {
      lwpointiterator_destroy(it);
      return LW_FALSE;
    }
  }
  lwpointiterator_destroy(it);
  return LW_TRUE;
}

Coverage Report

Created: 2026-08-13 07:11