[postgis-tickets] r15908 - Move simplify code to an in_place basis,

Paul Ramsey pramsey at cleverelephant.ca
Thu Oct 5 12:34:16 PDT 2017


Author: pramsey
Date: 2017-10-05 12:34:16 -0700 (Thu, 05 Oct 2017)
New Revision: 15908

Modified:
   trunk/liblwgeom/cunit/cu_algorithm.c
   trunk/liblwgeom/liblwgeom.h.in
   trunk/liblwgeom/liblwgeom_internal.h
   trunk/liblwgeom/lwcollection.c
   trunk/liblwgeom/lwgeom.c
   trunk/liblwgeom/lwgeom_api.c
   trunk/liblwgeom/lwline.c
   trunk/liblwgeom/lwpoly.c
   trunk/liblwgeom/ptarray.c
Log:
Move simplify code to an in_place basis, 
references #3877


Modified: trunk/liblwgeom/cunit/cu_algorithm.c
===================================================================
--- trunk/liblwgeom/cunit/cu_algorithm.c	2017-10-05 18:21:41 UTC (rev 15907)
+++ trunk/liblwgeom/cunit/cu_algorithm.c	2017-10-05 19:34:16 UTC (rev 15908)
@@ -927,59 +927,59 @@
 
 static void test_lwgeom_simplify(void)
 {
-		LWGEOM *l;
-		LWGEOM *g;
-		char *ewkt;
+	LWGEOM *l;
+	LWGEOM *g;
+	char *ewkt;
 
-		/* Simplify but only so far... */
-		g = lwgeom_from_wkt("LINESTRING(0 0, 1 0, 1 1, 0 1, 0 0)", LW_PARSER_CHECK_NONE);
-		l = lwgeom_simplify(g, 10, LW_TRUE);
-		ewkt = lwgeom_to_ewkt(l);
-		CU_ASSERT_STRING_EQUAL(ewkt, "LINESTRING(0 0,0 0)");
-		lwgeom_free(g);
-		lwgeom_free(l);
-		lwfree(ewkt);
+	/* Simplify but only so far... */
+	g = lwgeom_from_wkt("LINESTRING(0 0, 1 0, 1 1, 0 1, 0 0)", LW_PARSER_CHECK_NONE);
+	l = lwgeom_simplify(g, 10, LW_TRUE);
+	ewkt = lwgeom_to_ewkt(l);
+	CU_ASSERT_STRING_EQUAL(ewkt, "LINESTRING(0 0,0 0)");
+	lwgeom_free(g);
+	lwgeom_free(l);
+	lwfree(ewkt);
 
-		/* Simplify but only so far... */
-		g = lwgeom_from_wkt("POLYGON((0 0, 1 0, 1 1, 0 1, 0 0))", LW_PARSER_CHECK_NONE);
-		l = lwgeom_simplify(g, 10, LW_TRUE);
-		ewkt = lwgeom_to_ewkt(l);
-		CU_ASSERT_STRING_EQUAL(ewkt, "POLYGON((0 0,1 0,1 1,0 0))");
-		lwgeom_free(g);
-		lwgeom_free(l);
-		lwfree(ewkt);
+	/* Simplify but only so far... */
+	g = lwgeom_from_wkt("POLYGON((0 0, 1 0, 1 1, 0 1, 0 0))", LW_PARSER_CHECK_NONE);
+	l = lwgeom_simplify(g, 10, LW_TRUE);
+	ewkt = lwgeom_to_ewkt(l);
+	CU_ASSERT_STRING_EQUAL(ewkt, "POLYGON((0 0,1 0,1 1,0 0))");
+	lwgeom_free(g);
+	lwgeom_free(l);
+	lwfree(ewkt);
 
-		/* Simplify and collapse */
-		g = lwgeom_from_wkt("LINESTRING(0 0, 1 0, 1 1, 0 1, 0 0)", LW_PARSER_CHECK_NONE);
-		l = lwgeom_simplify(g, 10, LW_FALSE);
-		CU_ASSERT_EQUAL(l, NULL);
-		lwgeom_free(g);
-		lwgeom_free(l);
+	/* Simplify and collapse */
+	g = lwgeom_from_wkt("LINESTRING(0 0, 1 0, 1 1, 0 1, 0 0)", LW_PARSER_CHECK_NONE);
+	l = lwgeom_simplify(g, 10, LW_FALSE);
+	CU_ASSERT_EQUAL(l, NULL);
+	lwgeom_free(g);
+	lwgeom_free(l);
 
-		/* Simplify and collapse */
-		g = lwgeom_from_wkt("POLYGON((0 0, 1 0, 1 1, 0 1, 0 0))", LW_PARSER_CHECK_NONE);
-		l = lwgeom_simplify(g, 10, LW_FALSE);
-		CU_ASSERT_EQUAL(l, NULL);
-		lwgeom_free(g);
-		lwgeom_free(l);
+	/* Simplify and collapse */
+	g = lwgeom_from_wkt("POLYGON((0 0, 1 0, 1 1, 0 1, 0 0))", LW_PARSER_CHECK_NONE);
+	l = lwgeom_simplify(g, 10, LW_FALSE);
+	CU_ASSERT_EQUAL(l, NULL);
+	lwgeom_free(g);
+	lwgeom_free(l);
 
-		/* Not simplifiable */
-		g = lwgeom_from_wkt("LINESTRING(0 0, 50 1.00001, 100 0)", LW_PARSER_CHECK_NONE);
-		l = lwgeom_simplify(g, 1.0, LW_FALSE);
-		ewkt = lwgeom_to_ewkt(l);
-		CU_ASSERT_STRING_EQUAL(ewkt, "LINESTRING(0 0,50 1.00001,100 0)");
-		lwgeom_free(g);
-		lwgeom_free(l);
-		lwfree(ewkt);
+	/* Not simplifiable */
+	g = lwgeom_from_wkt("LINESTRING(0 0, 50 1.00001, 100 0)", LW_PARSER_CHECK_NONE);
+	l = lwgeom_simplify(g, 1.0, LW_FALSE);
+	ewkt = lwgeom_to_ewkt(l);
+	CU_ASSERT_STRING_EQUAL(ewkt, "LINESTRING(0 0,50 1.00001,100 0)");
+	lwgeom_free(g);
+	lwgeom_free(l);
+	lwfree(ewkt);
 
-		/* Simplifiable */
-		g = lwgeom_from_wkt("LINESTRING(0 0,50 0.99999,100 0)", LW_PARSER_CHECK_NONE);
-		l = lwgeom_simplify(g, 1.0, LW_FALSE);
-		ewkt = lwgeom_to_ewkt(l);
-		CU_ASSERT_STRING_EQUAL(ewkt, "LINESTRING(0 0,100 0)");
-		lwgeom_free(g);
-		lwgeom_free(l);
-		lwfree(ewkt);
+	/* Simplifiable */
+	g = lwgeom_from_wkt("LINESTRING(0 0,50 0.99999,100 0)", LW_PARSER_CHECK_NONE);
+	l = lwgeom_simplify(g, 1.0, LW_FALSE);
+	ewkt = lwgeom_to_ewkt(l);
+	CU_ASSERT_STRING_EQUAL(ewkt, "LINESTRING(0 0,100 0)");
+	lwgeom_free(g);
+	lwgeom_free(l);
+	lwfree(ewkt);
 }
 
 

Modified: trunk/liblwgeom/liblwgeom.h.in
===================================================================
--- trunk/liblwgeom/liblwgeom.h.in	2017-10-05 18:21:41 UTC (rev 15907)
+++ trunk/liblwgeom/liblwgeom.h.in	2017-10-05 19:34:16 UTC (rev 15908)
@@ -1224,6 +1224,8 @@
 
 
 
+extern LWGEOM* lwgeom_simplify(const LWGEOM *igeom, double dist, int preserve_collapsed);
+
 /****************************************************************
 * READ/WRITE FUNCTIONS
 *
@@ -1237,16 +1239,14 @@
 extern void lwgeom_reverse_in_place(LWGEOM *lwgeom);
 extern void lwgeom_force_clockwise(LWGEOM *lwgeom);
 extern void lwgeom_longitude_shift(LWGEOM *lwgeom);
-extern LWGEOM* lwgeom_simplify(const LWGEOM *igeom, double dist, int preserve_collapsed);
+extern void lwgeom_simplify_in_place(LWGEOM *igeom, double dist, int preserve_collapsed);
 extern void lwgeom_affine(LWGEOM *geom, const AFFINE *affine);
 extern void lwgeom_scale(LWGEOM *geom, const POINT4D *factors);
+extern LWGEOM* lwgeom_remove_repeated_points(const LWGEOM *in, double tolerance);
 
 
 
 
-/****************************************************************/
-
-
 /**
  * @brief wrap geometry on given cut x value
  *
@@ -1656,7 +1656,6 @@
 /**
 * Remove repeated points!
 */
-extern LWGEOM* lwgeom_remove_repeated_points(const LWGEOM *in, double tolerance);
 
 /**
  * Swap ordinate values in every vertex of the geometry.

Modified: trunk/liblwgeom/liblwgeom_internal.h
===================================================================
--- trunk/liblwgeom/liblwgeom_internal.h	2017-10-05 18:21:41 UTC (rev 15907)
+++ trunk/liblwgeom/liblwgeom_internal.h	2017-10-05 19:34:16 UTC (rev 15908)
@@ -201,10 +201,7 @@
 /**
  * @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, int preserve_collapsed);
-LWPOLY* lwpoly_simplify(const LWPOLY *ipoly, double dist, int preserve_collapsed);
-LWCOLLECTION* lwcollection_simplify(const LWCOLLECTION *igeom, double dist, int preserve_collapsed);
+void ptarray_simplify_in_place(POINTARRAY *pa, double epsilon, unsigned int minpts);
 
 /*
 * The possible ways a pair of segments can interact. Returned by lw_segment_intersects
@@ -375,6 +372,12 @@
 void ptarray_longitude_shift(POINTARRAY *pa);
 
 /*
+* Support for in place modification of point arrays, fast
+* function to move coordinate values around
+*/
+void ptarray_copy_point(POINTARRAY *pa, int from, int to);
+
+/*
 * Reverse
 */
 void ptarray_reverse_in_place(POINTARRAY *pa);

Modified: trunk/liblwgeom/lwcollection.c
===================================================================
--- trunk/liblwgeom/lwcollection.c	2017-10-05 18:21:41 UTC (rev 15907)
+++ trunk/liblwgeom/lwcollection.c	2017-10-05 19:34:16 UTC (rev 15908)
@@ -523,23 +523,7 @@
 	return v;
 }
 
-LWCOLLECTION* lwcollection_simplify(const LWCOLLECTION *igeom, double dist, int preserve_collapsed)
-{
- 	int i;
-	LWCOLLECTION *out = lwcollection_construct_empty(igeom->type, igeom->srid, FLAGS_GET_Z(igeom->flags), FLAGS_GET_M(igeom->flags));
 
-	if( lwcollection_is_empty(igeom) )
-		return out; /* should we return NULL instead ? */
-
-	for( i = 0; i < igeom->ngeoms; i++ )
-	{
-		LWGEOM *ngeom = lwgeom_simplify(igeom->geoms[i], dist, preserve_collapsed);
-		if ( ngeom ) out = lwcollection_add_lwgeom(out, ngeom);
-	}
-
-	return out;
-}
-
 int lwcollection_allows_subtype(int collectiontype, int subtype)
 {
 	if ( collectiontype == COLLECTIONTYPE )

Modified: trunk/liblwgeom/lwgeom.c
===================================================================
--- trunk/liblwgeom/lwgeom.c	2017-10-05 18:21:41 UTC (rev 15907)
+++ trunk/liblwgeom/lwgeom.c	2017-10-05 19:34:16 UTC (rev 15908)
@@ -1630,27 +1630,116 @@
 	}
 }
 
-LWGEOM* lwgeom_simplify(const LWGEOM *igeom, double dist, int preserve_collapsed)
+/**************************************************************/
+
+
+void
+lwgeom_simplify_in_place(LWGEOM *geom, double epsilon, int preserve_collapsed)
 {
-	switch (igeom->type)
+	switch (geom->type)
 	{
-	case POINTTYPE:
-	case MULTIPOINTTYPE:
-		return lwgeom_clone(igeom);
-	case LINETYPE:
-		return (LWGEOM*)lwline_simplify((LWLINE*)igeom, dist, preserve_collapsed);
-	case POLYGONTYPE:
-		return (LWGEOM*)lwpoly_simplify((LWPOLY*)igeom, dist, preserve_collapsed);
-	case MULTILINETYPE:
-	case MULTIPOLYGONTYPE:
-	case COLLECTIONTYPE:
-		return (LWGEOM*)lwcollection_simplify((LWCOLLECTION *)igeom, dist, preserve_collapsed);
-	default:
-		lwerror("%s: unsupported geometry type: %s", __func__, lwtype_name(igeom->type));
+		/* No-op! Cannot simplify points */
+		case POINTTYPE:
+			return;
+		case LINETYPE:
+		{
+			LWLINE *g = (LWLINE*)(geom);
+			POINTARRAY *pa = g->points;
+			ptarray_simplify_in_place(pa, epsilon, 2);
+			/* 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:
+		{
+			int 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;
+				ptarray_simplify_in_place(pa, epsilon, minpoints);
+				/* Drop collapsed rings */
+				if(pa->npoints < 4)
+				{
+					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 COLLECTIONTYPE:
+		{
+			int i, j = 0;
+			LWCOLLECTION *col = (LWCOLLECTION*)geom;
+			for (i = 0; i < col->ngeoms; i++)
+			{
+				LWGEOM *g = col->geoms[i];
+				if (!g) continue;
+				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;
+		}
 	}
-	return NULL;
+	return;
 }
 
+
+/**************************************************************/
+
+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)) return NULL;
+	return lwgeom_out;
+}
+
+
 double lwgeom_area(const LWGEOM *geom)
 {
 	int type = geom->type;

Modified: trunk/liblwgeom/lwgeom_api.c
===================================================================
--- trunk/liblwgeom/lwgeom_api.c	2017-10-05 18:21:41 UTC (rev 15907)
+++ trunk/liblwgeom/lwgeom_api.c	2017-10-05 19:34:16 UTC (rev 15908)
@@ -448,7 +448,43 @@
 	}
 }
 
+void
+ptarray_copy_point(POINTARRAY *pa, int from, int to)
+{
+	int ndims = FLAGS_NDIMS(pa->flags);
+	switch (ndims)
+	{
+		case 2:
+		{
+			POINT2D *p_from = (POINT2D*)(getPoint_internal(pa, from));
+			POINT2D *p_to = (POINT2D*)(getPoint_internal(pa, to));
+			*p_to = *p_from;
+			return;
+		}
+		case 3:
+		{
+			POINT3D *p_from = (POINT3D*)(getPoint_internal(pa, from));
+			POINT3D *p_to = (POINT3D*)(getPoint_internal(pa, to));
+			*p_to = *p_from;
+			return;
+		}
+		case 4:
+		{
+			POINT4D *p_from = (POINT4D*)(getPoint_internal(pa, from));
+			POINT4D *p_to = (POINT4D*)(getPoint_internal(pa, to));
+			*p_to = *p_from;
+			return;
+		}
+		default:
+		{
+			lwerror("%s: unsupported number of dimensions - %d", __func__, ndims);
+			return;
+		}
+	}
+	return;
+}
 
+
 /************************************************
  * debugging routines
  ************************************************/

Modified: trunk/liblwgeom/lwline.c
===================================================================
--- trunk/liblwgeom/lwline.c	2017-10-05 18:21:41 UTC (rev 15907)
+++ trunk/liblwgeom/lwline.c	2017-10-05 19:34:16 UTC (rev 15908)
@@ -532,44 +532,6 @@
 	return line->points->npoints;
 }
 
-LWLINE* lwline_simplify(const LWLINE *iline, double dist, int preserve_collapsed)
-{
-	static const int minvertices = 2; /* TODO: allow setting this */
-	LWLINE *oline;
-	POINTARRAY *pa;
-
-	LWDEBUG(2, "function called");
-
-	/* Skip empty case */
-	if( lwline_is_empty(iline) )
-		return NULL;
-
-	pa = ptarray_simplify(iline->points, dist, minvertices);
-	if ( ! pa ) return NULL;
-
-	/* Make sure single-point collapses have two points */
-	if ( pa->npoints == 1 )
-	{
-		/* Make sure single-point collapses have two points */
-		if ( preserve_collapsed )
-		{
-			POINT4D pt;
-			getPoint4d_p(pa, 0, &pt);
-			ptarray_append_point(pa, &pt, LW_TRUE);
-		}
-		/* Return null for collapse */
-		else
-		{
-			ptarray_free(pa);
-			return NULL;
-		}
-	}
-
-	oline = lwline_construct(iline->srid, NULL, pa);
-	oline->type = iline->type;
-	return oline;
-}
-
 double lwline_length(const LWLINE *line)
 {
 	if ( lwline_is_empty(line) )

Modified: trunk/liblwgeom/lwpoly.c
===================================================================
--- trunk/liblwgeom/lwpoly.c	2017-10-05 18:21:41 UTC (rev 15907)
+++ trunk/liblwgeom/lwpoly.c	2017-10-05 19:34:16 UTC (rev 15908)
@@ -450,62 +450,6 @@
 	return v;
 }
 
-LWPOLY* lwpoly_simplify(const LWPOLY *ipoly, double dist, int preserve_collapsed)
-{
-	int i;
-	LWPOLY *opoly = lwpoly_construct_empty(ipoly->srid, FLAGS_GET_Z(ipoly->flags), FLAGS_GET_M(ipoly->flags));
-
-	LWDEBUGF(2, "%s: simplifying polygon with %d rings", __func__, ipoly->nrings);
-
-	if ( lwpoly_is_empty(ipoly) )
-	{
-		lwpoly_free(opoly);
-		return NULL;
-	}
-
-	for ( i = 0; i < ipoly->nrings; i++ )
-	{
-		POINTARRAY *opts;
-		int minvertices = 0;
-
-		/* We'll still let holes collapse, but if we're preserving */
-		/* and this is a shell, we ensure it is kept */
-		if ( preserve_collapsed && i == 0 )
-			minvertices = 4;
-
-		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);
-
-		/* Less points than are needed to form a closed ring, we can't use this */
-		if ( opts->npoints < 4 )
-		{
-			LWDEBUGF(3, "ring%d skipped (% pts)", i, opts->npoints);
-			ptarray_free(opts);
-			if ( i ) continue;
-			else break; /* Don't scan holes if shell is collapsed */
-		}
-
-		/* Add ring to simplified polygon */
-		if( lwpoly_add_ring(opoly, opts) == LW_FAILURE )
-		{
-			lwpoly_free(opoly);
-			return NULL;
-		}
-	}
-
-	LWDEBUGF(3, "simplified polygon with %d rings", ipoly->nrings);
-	opoly->type = ipoly->type;
-
-	if( lwpoly_is_empty(opoly) )
-	{
-		lwpoly_free(opoly);
-		return NULL;
-	}
-
-	return opoly;
-}
-
 /**
 * Find the area of the outer ring - sum (area of inner rings).
 */

Modified: trunk/liblwgeom/ptarray.c
===================================================================
--- trunk/liblwgeom/ptarray.c	2017-10-05 18:21:41 UTC (rev 15907)
+++ trunk/liblwgeom/ptarray.c	2017-10-05 19:34:16 UTC (rev 15908)
@@ -1517,8 +1517,10 @@
 	return ptarray_remove_repeated_points_minpoints(in, tolerance, 2);
 }
 
+/************************************************************************/
+
 static void
-ptarray_dp_findsplit(POINTARRAY *pts, int p1, int p2, int *split, double *dist)
+ptarray_dp_findsplit_in_place(const POINTARRAY *pts, int p1, int p2, int *split, double *dist)
 {
 	int k;
 	const POINT2D *pk, *pa, *pb;
@@ -1556,75 +1558,105 @@
 			}
 		}
 		*dist = d;
-
-	} /* length---should be redone if can == 0 */
+	}
 	else
 	{
 		LWDEBUG(3, "segment too short, no split/no dist");
 		*dist = -1;
 	}
+}
 
+static int
+int_cmp(const void *a, const void *b)
+{
+	/* casting pointer types */
+    const int *ia = (const int *)a;
+    const int *ib = (const int *)b;
+	/* returns negative if b > a and positive if a > b */
+    return *ia - *ib;
 }
 
-POINTARRAY *
-ptarray_simplify(POINTARRAY *inpts, double epsilon, unsigned int minpts)
+void
+ptarray_simplify_in_place(POINTARRAY *pa, double epsilon, unsigned int minpts)
 {
-	int *stack;			/* recursion stack */
-	int sp=-1;			/* recursion stack pointer */
+	static size_t stack_size = 256;
+	int *stack, *outlist; /* recursion stack */
+	int stack_static[stack_size];
+	int outlist_static[stack_size];
+	int sp = -1; /* recursion stack pointer */
 	int p1, split;
+	int outn = 0;
+	int pai = 0;
+	int i;
 	double dist;
-	POINTARRAY *outpts;
-	POINT4D pt;
-
 	double eps_sqr = epsilon * epsilon;
 
-	/* Allocate recursion stack */
-	stack = lwalloc(sizeof(int)*inpts->npoints);
+	/* Do not try to simplify really short things */
+	if (pa->npoints < 3) return;
 
+	/* Only heap allocate book-keeping arrays if necessary */
+	if (pa->npoints > stack_size)
+	{
+		stack = lwalloc(sizeof(int) * pa->npoints);
+		outlist = lwalloc(sizeof(int) * pa->npoints);
+	}
+	else
+	{
+		stack = stack_static;
+		outlist = outlist_static;
+	}
+
 	p1 = 0;
-	stack[++sp] = inpts->npoints-1;
+	stack[++sp] = pa->npoints-1;
 
-	LWDEBUGF(2, "Input has %d pts and %d dims", inpts->npoints,
-	                                            FLAGS_NDIMS(inpts->flags));
-
-	/* Allocate output POINTARRAY, and add first point. */
-	outpts = ptarray_construct_empty(FLAGS_GET_Z(inpts->flags), FLAGS_GET_M(inpts->flags), inpts->npoints);
-	getPoint4d_p(inpts, 0, &pt);
-	ptarray_append_point(outpts, &pt, LW_FALSE);
-
-	LWDEBUG(3, "Added P0 to simplified point array (size 1)");
-
+	/* Add first point to output list */
+	outlist[outn++] = 0;
 	do
 	{
+		ptarray_dp_findsplit_in_place(pa, p1, stack[sp], &split, &dist);
 
-		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);
-
-		if (dist > eps_sqr || ( outpts->npoints+sp+1 < minpts && dist >= 0 ) )
+		if ((dist > eps_sqr) || ((outn + 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);
-			LWDEBUGF(4, "npoints , minpoints %d %d", outpts->npoints, minpts);
-			ptarray_append_point(outpts, &pt, LW_FALSE);
-
-			LWDEBUGF(4, "Added P%d to simplified point array (size: %d)", stack[sp], outpts->npoints);
-
+			outlist[outn++] = stack[sp];
 			p1 = stack[sp--];
 		}
+	}
+	while (!(sp<0));
 
-		LWDEBUGF(4, "stack pointer = %d", sp);
+	/* Put list of retained points into order */
+	qsort(outlist, outn, sizeof(int), int_cmp);
+	/* Copy retained points to front of array */
+	for (i = 0; i < outn; i++)
+	{
+		int j = outlist[i];
+		/* Indexes the same, means no copy required */
+		if (j == pai)
+		{
+			pai++;
+			continue;
+		}
+		/* Indexes different, copy value down */
+		ptarray_copy_point(pa, j, pai++);
 	}
-	while (! (sp<0) );
 
-	lwfree(stack);
-	return outpts;
+	/* Adjust point count on array */
+	pa->npoints = outn;
+
+	/* Only free if arrays are on heap */
+	if (stack != stack_static)
+		lwfree(stack);
+	if (outlist != outlist_static)
+		lwfree(outlist);
+
+	return;
 }
 
+/************************************************************************/
+
 /**
 * Find the 2d length of the given #POINTARRAY, using circular
 * arc interpolation between each coordinate triple.



More information about the postgis-tickets mailing list