added st_simplifypreservestructure

It allows to avoid removing features because of the simplification

SELECT '9', ST_astext(ST_Simplifyextended('POLYGON((0 0, 0 0,  0 0, 0 0))', 10, 1));

 ?column? |         st_astext
----------+----------------------------
 9        | POLYGON((0 0,0 0,0 0,0 0))

liblwgeom/liblwgeom.h.in

@@ -927,7 +927,7 @@ extern LWGEOM* lwgeom_force_3dz(const LWGEOM *geom);
 extern LWGEOM* lwgeom_force_3dm(const LWGEOM *geom);
 extern LWGEOM* lwgeom_force_4d(const LWGEOM *geom);
 
-extern LWGEOM* lwgeom_simplify(const LWGEOM *igeom, double dist);
+extern LWGEOM* lwgeom_simplify(const LWGEOM *igeom, double dist, uint8_t flags);
 extern LWGEOM* lwgeom_set_effective_area(const LWGEOM *igeom, double area);
 
 /* 
@@ -1842,6 +1842,9 @@ LWGEOM_UNPARSER_RESULT;
 #define TWKB_NO_TYPE 0x08 /*No type because it is a sub geoemtry*/
 #define TWKB_NO_ID 0x10 /*No ID because it is a subgeoemtry*/
 
+/* simplification flags */
+#define SIM_PRESERVE_GEOM 0x01
+
 /*
 ** New parsing and unparsing functions.
 */

liblwgeom/liblwgeom_internal.h

@@ -185,9 +185,9 @@ extern int32_t lw_get_int32_t(const uint8_t *loc);
  * @param minpts minimun number of points to retain, if possible.
  */
 POINTARRAY* ptarray_simplify(POINTARRAY *inpts, double epsilon, unsigned int minpts);
-LWLINE* lwline_simplify(const LWLINE *iline, double dist);
-LWPOLY* lwpoly_simplify(const LWPOLY *ipoly, double dist);
-LWCOLLECTION* lwcollection_simplify(const LWCOLLECTION *igeom, double dist);
+LWLINE* lwline_simplify(const LWLINE *iline, double dist, uint8_t flags);
+LWPOLY* lwpoly_simplify(const LWPOLY *ipoly, double dist, uint8_t flags);
+LWCOLLECTION* lwcollection_simplify(const LWCOLLECTION *igeom, double dist, uint8_t flags);
 
 /*
 * Computational geometry

liblwgeom/lwcollection.c

@@ -503,7 +503,7 @@ int lwcollection_count_vertices(LWCOLLECTION *col)
 	return v;
 }
 
-LWCOLLECTION* lwcollection_simplify(const LWCOLLECTION *igeom, double dist)
+LWCOLLECTION* lwcollection_simplify(const LWCOLLECTION *igeom, double dist, uint8_t flags)
 {
  	int i;
 	LWCOLLECTION *out = lwcollection_construct_empty(igeom->type, igeom->srid, FLAGS_GET_Z(igeom->flags), FLAGS_GET_M(igeom->flags));
@@ -513,7 +513,7 @@ LWCOLLECTION* lwcollection_simplify(const LWCOLLECTION *igeom, double dist)
 
 	for( i = 0; i < igeom->ngeoms; i++ )
 	{
-		LWGEOM *ngeom = lwgeom_simplify(igeom->geoms[i], dist);
+		LWGEOM *ngeom = lwgeom_simplify(igeom->geoms[i], dist, flags);
 		if ( ngeom ) out = lwcollection_add_lwgeom(out, ngeom);
 	}
 

liblwgeom/lwgeom.c

@@ -1511,21 +1511,21 @@ void lwgeom_set_srid(LWGEOM *geom, int32_t srid)
 	}
 }
 
-LWGEOM* lwgeom_simplify(const LWGEOM *igeom, double dist)
+LWGEOM* lwgeom_simplify(const LWGEOM *igeom, double dist, uint8_t flags) 
 {
 	switch (igeom->type)
 	{
 	case POINTTYPE:
 	case MULTIPOINTTYPE:
 		return lwgeom_clone(igeom);
 	case LINETYPE:
-		return (LWGEOM*)lwline_simplify((LWLINE*)igeom, dist);
+		return (LWGEOM*)lwline_simplify((LWLINE*)igeom, dist, flags);
 	case POLYGONTYPE:
-		return (LWGEOM*)lwpoly_simplify((LWPOLY*)igeom, dist);
+		return (LWGEOM*)lwpoly_simplify((LWPOLY*)igeom, dist, flags);
 	case MULTILINETYPE:
 	case MULTIPOLYGONTYPE:
 	case COLLECTIONTYPE:
-		return (LWGEOM*)lwcollection_simplify((LWCOLLECTION *)igeom, dist);
+		return (LWGEOM*)lwcollection_simplify((LWCOLLECTION *)igeom, dist, flags);
 	default:
 		lwerror("lwgeom_simplify: unsupported geometry type: %s",lwtype_name(igeom->type));
 	}

liblwgeom/lwline.c

@@ -475,7 +475,7 @@ int lwline_count_vertices(LWLINE *line)
 	return line->points->npoints;
 }
 
-LWLINE* lwline_simplify(const LWLINE *iline, double dist)
+LWLINE* lwline_simplify(const LWLINE *iline, double dist, uint8_t flags)
 {
 	LWLINE *oline;
 
@@ -485,7 +485,10 @@ LWLINE* lwline_simplify(const LWLINE *iline, double dist)
 	if( lwline_is_empty(iline) )
 		return lwline_clone(iline);
 
-	static const int minvertices = 0; /* TODO: allow setting this */
+	int minvertices = 0; /* TODO: allow setting this */
+	if (flags & SIM_PRESERVE_GEOM) {
+		minvertices = 2;
+	}
 	oline = lwline_construct(iline->srid, NULL, ptarray_simplify(iline->points, dist, minvertices));
 	oline->type = iline->type;
 	return oline;

liblwgeom/lwpoly.c

@@ -348,9 +348,16 @@ int lwpoly_count_vertices(LWPOLY *poly)
 	return v;
 }
 
-LWPOLY* lwpoly_simplify(const LWPOLY *ipoly, double dist)
+LWPOLY* lwpoly_simplify(const LWPOLY *ipoly, double dist, uint8_t flags)
 {
 	int i;
+	int minvertices = 0;
+
+	if (flags & SIM_PRESERVE_GEOM) {
+		minvertices = 4;
+		LWDEBUG(2, "preserving geometry");
+	}
+
 	LWPOLY *opoly = lwpoly_construct_empty(ipoly->srid, FLAGS_GET_Z(ipoly->flags), FLAGS_GET_M(ipoly->flags));
 
 	LWDEBUGF(2, "simplify_polygon3d: simplifying polygon with %d rings", ipoly->nrings);
@@ -360,7 +367,6 @@ LWPOLY* lwpoly_simplify(const LWPOLY *ipoly, double dist)
 
 	for (i = 0; i < ipoly->nrings; i++)
 	{
-		static const int minvertices = 0; /* TODO: allow setting this */
 		POINTARRAY *opts = ptarray_simplify(ipoly->rings[i], dist, minvertices);
 
 		LWDEBUGF(3, "ring%d simplified from %d to %d points", i, ipoly->rings[i]->npoints, opts->npoints);

liblwgeom/ptarray.c

@@ -1482,7 +1482,7 @@ POINTARRAY *
 ptarray_simplify(POINTARRAY *inpts, double epsilon, unsigned int minpts)
 {
 	int *stack;			/* recursion stack */
-	int sp=-1;			/* recursion stack pointer */
+	int sp=0;			/* recursion stack pointer */
 	int p1, split;
 	double dist;
 	POINTARRAY *outpts;
@@ -1491,11 +1491,15 @@ ptarray_simplify(POINTARRAY *inpts, double epsilon, unsigned int minpts)
 	/* Allocate recursion stack */
 	stack = lwalloc(sizeof(int)*inpts->npoints);
 
+	if (inpts->npoints < minpts) {
+		lwerror("min number of points is greater than geometry points");
+	}
+
 	p1 = 0;
-	stack[++sp] = inpts->npoints-1;
+	stack[sp] = inpts->npoints-1;
 
-	LWDEBUGF(2, "Input has %d pts and %d dims", inpts->npoints,
-	                                            FLAGS_NDIMS(inpts->flags));
+	LWDEBUGF(2, "Input has %d pts and %d dims. Min poinds %d", inpts->npoints,
+	                                            FLAGS_NDIMS(inpts->flags), minpts);
 
 	/* Allocate output POINTARRAY, and add first point. */
 	outpts = ptarray_construct_empty(FLAGS_GET_Z(inpts->flags), FLAGS_GET_M(inpts->flags), inpts->npoints);
@@ -1509,17 +1513,18 @@ ptarray_simplify(POINTARRAY *inpts, double epsilon, unsigned int minpts)
 
 		ptarray_dp_findsplit(inpts, p1, stack[sp], &split, &dist);
 
-		LWDEBUGF(3, "Farthest point from P%d-P%d is P%d (dist. %g)", p1, stack[sp], split, dist);
+		LWDEBUGF(3, "Farthest point from P%d-P%d is P%d (dist. %g) (expected output: %d, %d)", p1, stack[sp], split, dist, outpts->npoints+sp +1, minpts);
 
-		if (dist > epsilon || ( outpts->npoints+sp+1 < minpts && dist > 0 ) )
+		// dist == -1 when there are only two points
+		if (dist > epsilon || ( outpts->npoints+sp+1 < minpts && dist >= 0 ) )
 		{
 			LWDEBUGF(4, "Added P%d to stack (outpts:%d)", split, sp);
 			stack[++sp] = split;
 		}
 		else
 		{
 			getPoint4d_p(inpts, stack[sp], &pt);
-			ptarray_append_point(outpts, &pt, LW_FALSE);
+			ptarray_append_point(outpts, &pt, outpts->npoints < minpts ? LW_TRUE : LW_FALSE);
 
 			LWDEBUGF(4, "Added P%d to simplified point array (size: %d)", stack[sp], outpts->npoints);
 

postgis/lwgeom_functions_analytic.c

@@ -31,6 +31,7 @@
 
 /* Prototypes */
 Datum LWGEOM_simplify2d(PG_FUNCTION_ARGS);
+Datum LWGEOM_SimplifyPreserveStructure(PG_FUNCTION_ARGS);
 Datum LWGEOM_SetEffectiveArea(PG_FUNCTION_ARGS);
 Datum ST_LineCrossingDirection(PG_FUNCTION_ARGS);
 
@@ -56,7 +57,38 @@ Datum LWGEOM_simplify2d(PG_FUNCTION_ARGS)
 	dist = PG_GETARG_FLOAT8(1);
 	in = lwgeom_from_gserialized(geom);
 
-	out = lwgeom_simplify(in, dist);
+	out = lwgeom_simplify(in, dist, 0);
+	if ( ! out ) PG_RETURN_NULL();
+
+	/* COMPUTE_BBOX TAINTING */
+	if ( in->bbox ) lwgeom_add_bbox(out);
+
+	result = geometry_serialize(out);
+	lwgeom_free(out);
+	PG_FREE_IF_COPY(geom, 0);
+	PG_RETURN_POINTER(result);
+}
+
+PG_FUNCTION_INFO_V1(LWGEOM_SimplifyPreserveStructure);
+Datum LWGEOM_SimplifyPreserveStructure(PG_FUNCTION_ARGS)
+{
+	GSERIALIZED *geom = (GSERIALIZED *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
+	GSERIALIZED *result;
+	int type = gserialized_get_type(geom);
+	LWGEOM *in;
+	LWGEOM *out;
+	double dist;
+	uint8_t flags;
+
+  if ( type == POINTTYPE || type == MULTIPOINTTYPE )
+    PG_RETURN_POINTER(geom);
+
+	dist = PG_GETARG_FLOAT8(1);
+	in = lwgeom_from_gserialized(geom);
+
+	dist = PG_GETARG_DATUM(1);
+	flags = PG_GETARG_INT32(2);
+	out = lwgeom_simplify(in, dist, flags);
 	if ( ! out ) PG_RETURN_NULL();
 
 	/* COMPUTE_BBOX TAINTING */
@@ -73,7 +105,7 @@ Datum LWGEOM_SetEffectiveArea(PG_FUNCTION_ARGS)
 {
 	GSERIALIZED *geom = (GSERIALIZED *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
 	GSERIALIZED *result;
-  int type = gserialized_get_type(geom);
+	int type = gserialized_get_type(geom);
 	LWGEOM *in;
 	LWGEOM *out;
 	double area;

postgis/postgis.sql.in

@@ -2844,6 +2844,12 @@ CREATE OR REPLACE FUNCTION ST_Simplify(geometry, float8)
 	RETURNS geometry
 	AS 'MODULE_PATHNAME', 'LWGEOM_simplify2d'
 	LANGUAGE 'c' IMMUTABLE STRICT;
+
+-- Availability: 2.2.0
+CREATE OR REPLACE FUNCTION ST_SimplifyPreserveStructure(geometry, float8, flags int4 DEFAULT 1)
+	RETURNS geometry
+	AS 'MODULE_PATHNAME', 'LWGEOM_SimplifyExtended2d'
+	LANGUAGE 'c' IMMUTABLE STRICT;
 
 -- Availability: 2.2.0
 CREATE OR REPLACE FUNCTION ST_SetEffectiveArea(geometry,  float8 default 0)