[postgis-tickets] r15925 - Move remove_repeated_points to an in_place basis

Paul Ramsey pramsey at cleverelephant.ca
Fri Oct 6 14:19:45 PDT 2017


Author: pramsey
Date: 2017-10-06 14:19:44 -0700 (Fri, 06 Oct 2017)
New Revision: 15925

Modified:
   trunk/liblwgeom/cunit/cu_algorithm.c
   trunk/liblwgeom/liblwgeom_internal.h
   trunk/liblwgeom/lwcollection.c
   trunk/liblwgeom/lwgeom.c
   trunk/liblwgeom/lwline.c
   trunk/liblwgeom/lwmpoint.c
   trunk/liblwgeom/lwpoly.c
   trunk/liblwgeom/ptarray.c
Log:
Move remove_repeated_points to an in_place basis
References #3877


Modified: trunk/liblwgeom/cunit/cu_algorithm.c
===================================================================
--- trunk/liblwgeom/cunit/cu_algorithm.c	2017-10-06 13:31:58 UTC (rev 15924)
+++ trunk/liblwgeom/cunit/cu_algorithm.c	2017-10-06 21:19:44 UTC (rev 15925)
@@ -1003,6 +1003,28 @@
 	CU_ASSERT_EQUAL(gh, rs);
 }
 
+static void test_lwgeom_remote_repeated_points(void)
+{
+	LWGEOM *g;
+	char *ewkt;
+
+	// g = lwgeom_from_wkt("MULTIPOINT(0 0, 10 0, 10 10, 10 10, 0 10, 0 10, 0 10, 0 0, 0 0, 0 0, 5 5, 0 0, 5 5)", LW_PARSER_CHECK_NONE);
+	// lwgeom_remove_repeated_points_in_place(g, 1);
+	// ewkt = lwgeom_to_ewkt(g);
+	// CU_ASSERT_STRING_EQUAL(ewkt, "MULTIPOINT(0 0,10 0,10 10,0 10,5 5)");
+	// printf("%s\n", ewkt);
+	// lwgeom_free(g);
+	// lwfree(ewkt);
+
+	g = lwgeom_from_wkt("MULTIPOINT(0 0,10 0,10 10,10 10,0 10,0 10,0 10,0 0,0 0,0 0,5 5,5 5,5 8,8 8,8 8,8 8,8 5,8 5,5 5,5 5,5 5,5 5,5 5,50 50,50 50,50 50,50 60,50 60,50 60,60 60,60 50,60 50,50 50,55 55,55 58,58 58,58 55,58 55,55 55)", LW_PARSER_CHECK_NONE);
+	lwgeom_remove_repeated_points_in_place(g, 1);
+	ewkt = lwgeom_to_ewkt(g);
+	CU_ASSERT_STRING_EQUAL(ewkt, "MULTIPOINT(0 0,10 0,10 10,0 10,5 5,5 8,8 8,8 5,50 50,50 60,60 60,60 50,55 55,55 58,58 58,58 55)");
+	printf("%s\n", ewkt);
+	lwgeom_free(g);
+	lwfree(ewkt);
+}
+
 static void test_lwgeom_simplify(void)
 {
 	LWGEOM *l;
@@ -1292,4 +1314,5 @@
 	PG_ADD_TEST(suite,test_median_handles_3d_correctly);
 	PG_ADD_TEST(suite,test_median_robustness);
 	PG_ADD_TEST(suite,test_lwpoly_construct_circle);
+	PG_ADD_TEST(suite,test_lwgeom_remote_repeated_points);
 }

Modified: trunk/liblwgeom/liblwgeom_internal.h
===================================================================
--- trunk/liblwgeom/liblwgeom_internal.h	2017-10-06 13:31:58 UTC (rev 15924)
+++ trunk/liblwgeom/liblwgeom_internal.h	2017-10-06 21:19:44 UTC (rev 15925)
@@ -398,12 +398,10 @@
 /*
 * Repeated points
 */
-POINTARRAY *ptarray_remove_repeated_points_minpoints(const POINTARRAY *in, double tolerance, int minpoints);
 POINTARRAY *ptarray_remove_repeated_points(const POINTARRAY *in, double tolerance);
-LWGEOM* lwmpoint_remove_repeated_points(const LWMPOINT *in, double tolerance);
 LWGEOM* lwline_remove_repeated_points(const LWLINE *in, double tolerance);
-LWGEOM* lwcollection_remove_repeated_points(const LWCOLLECTION *in, double tolerance);
-LWGEOM* lwpoly_remove_repeated_points(const LWPOLY *in, double tolerance);
+void ptarray_remove_repeated_points_in_place(POINTARRAY *pa, double tolerance, int min_points);
+void lwgeom_remove_repeated_points_in_place(LWGEOM *geom, double tolerance);
 
 /*
 * Closure test

Modified: trunk/liblwgeom/lwcollection.c
===================================================================
--- trunk/liblwgeom/lwcollection.c	2017-10-06 13:31:58 UTC (rev 15924)
+++ trunk/liblwgeom/lwcollection.c	2017-10-06 21:19:44 UTC (rev 15925)
@@ -456,24 +456,7 @@
 	return outcol;
 }
 
-LWGEOM*
-lwcollection_remove_repeated_points(const LWCOLLECTION *coll, double tolerance)
-{
-	uint32_t i;
-	LWGEOM **newgeoms;
 
-	newgeoms = lwalloc(sizeof(LWGEOM *)*coll->ngeoms);
-	for (i=0; i<coll->ngeoms; i++)
-	{
-		newgeoms[i] = lwgeom_remove_repeated_points(coll->geoms[i], tolerance);
-	}
-
-	return (LWGEOM*)lwcollection_construct(coll->type,
-	                                       coll->srid, coll->bbox ? gbox_copy(coll->bbox) : NULL,
-	                                       coll->ngeoms, newgeoms);
-}
-
-
 LWCOLLECTION*
 lwcollection_force_dims(const LWCOLLECTION *col, int hasz, int hasm)
 {

Modified: trunk/liblwgeom/lwgeom.c
===================================================================
--- trunk/liblwgeom/lwgeom.c	2017-10-06 13:31:58 UTC (rev 15924)
+++ trunk/liblwgeom/lwgeom.c	2017-10-06 21:19:44 UTC (rev 15925)
@@ -1486,53 +1486,9 @@
 
 extern LWGEOM* lwgeom_remove_repeated_points(const LWGEOM *in, double tolerance)
 {
-	LWDEBUGF(4, "lwgeom_remove_repeated_points got type %s",
-	         lwtype_name(in->type));
-
-	if(lwgeom_is_empty(in))
-	{
-		return lwgeom_clone_deep(in);
-	}
-
-	switch (in->type)
-	{
-	case MULTIPOINTTYPE:
-		return lwmpoint_remove_repeated_points((LWMPOINT*)in, tolerance);
-		break;
-	case LINETYPE:
-		return lwline_remove_repeated_points((LWLINE*)in, tolerance);
-
-	case MULTILINETYPE:
-	case COLLECTIONTYPE:
-	case MULTIPOLYGONTYPE:
-	case POLYHEDRALSURFACETYPE:
-		return lwcollection_remove_repeated_points((LWCOLLECTION *)in, tolerance);
-
-	case POLYGONTYPE:
-		return lwpoly_remove_repeated_points((LWPOLY *)in, tolerance);
-		break;
-
-	case POINTTYPE:
-	case TRIANGLETYPE:
-	case TINTYPE:
-		/* No point is repeated for a single point, or for Triangle or TIN */
-		return lwgeom_clone_deep(in);
-
-	case CIRCSTRINGTYPE:
-	case COMPOUNDTYPE:
-	case MULTICURVETYPE:
-	case CURVEPOLYTYPE:
-	case MULTISURFACETYPE:
-		/* Dunno how to handle these, will return untouched */
-		return lwgeom_clone_deep(in);
-
-	default:
-		lwnotice("%s: unsupported geometry type: %s",
-		         __func__, lwtype_name(in->type));
-		return lwgeom_clone_deep(in);
-		break;
-	}
-	return 0;
+	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)
@@ -1624,10 +1580,152 @@
 	}
 }
 
+
 /**************************************************************/
 
 
 void
+lwgeom_remove_repeated_points_in_place(LWGEOM *geom, double tolerance)
+{
+	switch (geom->type)
+	{
+		/* No-op! Cannot remote points */
+		case POINTTYPE:
+			return;
+		case LINETYPE:
+		{
+			LWLINE *g = (LWLINE*)(geom);
+			POINTARRAY *pa = g->points;
+			ptarray_remove_repeated_points_in_place(pa, tolerance, 2);
+			/* Invalid output */
+			if (pa->npoints == 1 && pa->maxpoints > 1)
+			{
+				/* Use first point as last point */
+				pa->npoints = 2;
+				ptarray_copy_point(pa, 0, 1);
+			}
+			break;
+		}
+		case POLYGONTYPE:
+		{
+			int i, j = 0;
+			LWPOLY *g = (LWPOLY*)(geom);
+			for (i = 0; i < g->nrings; i++)
+			{
+				POINTARRAY *pa = g->rings[i];
+				int minpoints = 4;
+				/* Skip zero'ed out rings */
+				if(!pa)
+					continue;
+				ptarray_remove_repeated_points_in_place(pa, tolerance, minpoints);
+				/* Drop collapsed rings */
+				if(pa->npoints < 4)
+				{
+					ptarray_free(pa);
+					continue;
+				}
+				g->rings[j++] = pa;
+			}
+			/* Update ring count */
+			g->nrings = j;
+			break;
+		}
+		case MULTIPOINTTYPE:
+		{
+			static int out_stack_size = 32;
+			double tolsq = tolerance*tolerance;
+			int i, j, n = 0;
+			LWMPOINT *mpt = (LWMPOINT *)(geom);
+			LWPOINT **out;
+			LWPOINT *out_stack[out_stack_size];
+			int use_heap = (mpt->ngeoms > out_stack_size);
+
+			/* No-op on empty */
+			if (mpt->ngeoms == 0) return;
+
+			/* We cannot write directly back to the multipoint */
+			/* geoms array because we're reading out of it still */
+			/* so we use a side array */
+			if (use_heap)
+				out = lwalloc(sizeof(LWMPOINT *) * mpt->ngeoms);
+			else
+				out = out_stack;
+
+			/* Inefficient O(n^2) implementation */
+			for (i = 0; i < mpt->ngeoms; i++)
+			{
+				int seen = 0;
+				LWPOINT *p1 = mpt->geoms[i];
+				const POINT2D *pt1 = getPoint2d_cp(p1->point, 0);
+				for (j = 0; j < n; j++)
+				{
+					LWPOINT *p2 = out[j];
+					const POINT2D *pt2 = getPoint2d_cp(p2->point, 0);
+					if (distance2d_sqr_pt_pt(pt1, pt2) <= tolsq)
+					{
+						seen = 1;
+						break;
+					}
+				}
+				if (seen) continue;
+				out[n++] = p1;
+			}
+
+			/* Copy remaining points back into the input */
+			/* array */
+			memcpy(mpt->geoms, out, sizeof(LWPOINT *) * n);
+			mpt->ngeoms = n;
+			if (use_heap) lwfree(out);
+			return;
+		}
+
+		case CIRCSTRINGTYPE:
+			/* Dunno how to handle these, will return untouched */
+			return;
+
+		/* Can process most multi* types as generic collection */
+		case MULTILINETYPE:
+		case MULTIPOLYGONTYPE:
+		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:
+		{
+			int i, j = 0;
+			LWCOLLECTION *col = (LWCOLLECTION*)(geom);
+			for (i = 0; i < col->ngeoms; i++)
+			{
+				LWGEOM *g = col->geoms[i];
+				if (!g) continue;
+				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;
+		}
+	}
+	return;
+}
+
+
+/**************************************************************/
+
+void
 lwgeom_simplify_in_place(LWGEOM *geom, double epsilon, int preserve_collapsed)
 {
 	switch (geom->type)
@@ -1722,9 +1820,6 @@
 	return;
 }
 
-
-/**************************************************************/
-
 LWGEOM* lwgeom_simplify(const LWGEOM *igeom, double dist, int preserve_collapsed)
 {
 	LWGEOM *lwgeom_out = lwgeom_clone_deep(igeom);
@@ -1733,7 +1828,9 @@
 	return lwgeom_out;
 }
 
+/**************************************************************/
 
+
 double lwgeom_area(const LWGEOM *geom)
 {
 	int type = geom->type;

Modified: trunk/liblwgeom/lwline.c
===================================================================
--- trunk/liblwgeom/lwline.c	2017-10-06 13:31:58 UTC (rev 15924)
+++ trunk/liblwgeom/lwline.c	2017-10-06 21:19:44 UTC (rev 15925)
@@ -449,13 +449,7 @@
 LWGEOM*
 lwline_remove_repeated_points(const LWLINE *lwline, double tolerance)
 {
-	POINTARRAY* npts = ptarray_remove_repeated_points_minpoints(lwline->points, tolerance, 2);
-
-	LWDEBUGF(3, "%s: npts %p", __func__, npts);
-
-	return (LWGEOM*)lwline_construct(lwline->srid,
-	                                 lwline->bbox ? gbox_copy(lwline->bbox) : 0,
-	                                 npts);
+	return lwgeom_remove_repeated_points((LWGEOM*)lwline, tolerance);
 }
 
 int

Modified: trunk/liblwgeom/lwmpoint.c
===================================================================
--- trunk/liblwgeom/lwmpoint.c	2017-10-06 13:31:58 UTC (rev 15924)
+++ trunk/liblwgeom/lwmpoint.c	2017-10-06 21:19:44 UTC (rev 15925)
@@ -88,41 +88,7 @@
 	lwfree(mpt);
 }
 
-LWGEOM*
-lwmpoint_remove_repeated_points(const LWMPOINT *mpoint, double tolerance)
-{
-	uint32_t nnewgeoms;
-	uint32_t i, j;
-	LWGEOM **newgeoms;
-	LWGEOM *lwpt1, *lwpt2;
 
-	newgeoms = lwalloc(sizeof(LWGEOM *)*mpoint->ngeoms);
-	nnewgeoms = 0;
-	for (i=0; i<mpoint->ngeoms; ++i)
-	{
-		lwpt1 = (LWGEOM*)mpoint->geoms[i];
-		/* Brute force, may be optimized by building an index */
-		int seen=0;
-		for (j=0; j<nnewgeoms; ++j)
-		{
-			lwpt2 = (LWGEOM*)newgeoms[j];
-			if ( lwgeom_mindistance2d(lwpt1, lwpt2) <= tolerance )
-			{
-				seen=1;
-				break;
-			}
-		}
-		if ( seen ) continue;
-		newgeoms[nnewgeoms++] = lwgeom_clone_deep(lwpt1);
-	}
-
-	return (LWGEOM*)lwcollection_construct(mpoint->type,
-	                                       mpoint->srid,
-										   mpoint->bbox ? gbox_copy(mpoint->bbox) : NULL,
-	                                       nnewgeoms, newgeoms);
-
-}
-
 LWMPOINT*
 lwmpoint_from_lwgeom(const LWGEOM *g)
 {

Modified: trunk/liblwgeom/lwpoly.c
===================================================================
--- trunk/liblwgeom/lwpoly.c	2017-10-06 13:31:58 UTC (rev 15924)
+++ trunk/liblwgeom/lwpoly.c	2017-10-06 21:19:44 UTC (rev 15925)
@@ -387,25 +387,6 @@
 	return ret;
 }
 
-LWGEOM*
-lwpoly_remove_repeated_points(const LWPOLY *poly, double tolerance)
-{
-	uint32_t i;
-	POINTARRAY **newrings;
-
-	newrings = lwalloc(sizeof(POINTARRAY *)*poly->nrings);
-	for (i=0; i<poly->nrings; i++)
-	{
-		newrings[i] = ptarray_remove_repeated_points_minpoints(poly->rings[i], tolerance, 4);
-	}
-
-	return (LWGEOM*)lwpoly_construct(poly->srid,
-	                                 poly->bbox ? gbox_copy(poly->bbox) : NULL,
-	                                 poly->nrings, newrings);
-
-}
-
-
 LWPOLY*
 lwpoly_force_dims(const LWPOLY *poly, int hasz, int hasm)
 {

Modified: trunk/liblwgeom/ptarray.c
===================================================================
--- trunk/liblwgeom/ptarray.c	2017-10-06 13:31:58 UTC (rev 15924)
+++ trunk/liblwgeom/ptarray.c	2017-10-06 21:19:44 UTC (rev 15925)
@@ -1429,93 +1429,86 @@
  *
  * Always returns a newly allocated object.
  */
-POINTARRAY *
+static POINTARRAY *
 ptarray_remove_repeated_points_minpoints(const POINTARRAY *in, double tolerance, int minpoints)
 {
-	POINTARRAY *out;
-	size_t ptsize = ptarray_point_size(in);
-	int has_z = FLAGS_GET_Z(in->flags);
-	int has_m = FLAGS_GET_M(in->flags);
-	int ipn = 1, opn = 1;
-	const POINT2D *prev_point, *this_point;
-	uint8_t *p1, *p2;
-	double tolsq = tolerance * tolerance;
+	POINTARRAY *out = ptarray_clone_deep(in);
+	ptarray_remove_repeated_points_in_place(out, tolerance, minpoints);
+	return out;
+}
 
-	LWDEBUGF(3, "%s called", __func__);
+POINTARRAY *
+ptarray_remove_repeated_points(const POINTARRAY *in, double tolerance)
+{
+	return ptarray_remove_repeated_points_minpoints(in, tolerance, 2);
+}
 
-	/* Single or zero point arrays can't have duplicates */
-	if ( in->npoints < 3 ) return ptarray_clone_deep(in);
 
-	/* Condition minpoints */
-	if ( minpoints < 1 ) minpoints = 1;
+void
+ptarray_remove_repeated_points_in_place(POINTARRAY *pa, double tolerance, int min_points)
+{
+	int i;
+	double tolsq = tolerance * tolerance;
+	const POINT2D *last = NULL;
+	const POINT2D *pt;
+	int n_points = pa->npoints;
+	int n_points_out = 0;
+	int pt_size = ptarray_point_size(pa);
+	double dsq;
 
-	/* Allocate enough output space for all points */
-	out = ptarray_construct(has_z, has_m, in->npoints);
+	/* No-op on short inputs */
+	if ( n_points <= 2 ) return;
 
-	/* Keep the first point */
-	p1 = getPoint_internal(out, 0);
-	p2 = getPoint_internal(in, 0);
-	memcpy(p1, p2, ptsize);
-
-	/* Now fill up the actual points */
-	prev_point = getPoint2d_cp(in, 0);
-	LWDEBUGF(3, " first point copied, out points: %d", opn);
-	for ( ipn = 1; ipn < in->npoints; ipn++ )
+	for (i = 0; i < n_points; i++)
 	{
-		this_point = getPoint2d_cp(in, ipn);
+		int last_point = (i == n_points-1);
 
-		/*
-		* If number of points left > number of points we need
-		* then it's still OK to drop dupes
-		*/
-		if ( in->npoints - ipn > minpoints - opn )
+		/* Look straight into the abyss */
+		pt = getPoint2d_cp(pa, i);
+
+		/* Preserve first point always */
+		if (last)
 		{
-			if (tolerance > 0.0)
+			/* Don't drop points if we are running short of points */
+	        if (n_points - i > min_points - n_points_out)
 			{
-				/* within the removal tolerance? */
-				double dsq = distance2d_sqr_pt_pt(prev_point, this_point);
-				if (dsq <= tolsq)
-					/* then skip it */
-					continue;
+				if (tolerance > 0.0)
+				{
+					/* Only drop points that are within our tolerance */
+					dsq = distance2d_sqr_pt_pt(last, pt);
+					/* Allow any point but the last one to be dropped */
+					if (!last_point && dsq <= tolsq)
+					{
+						continue;
+					}
+				}
+				else
+				{
+					/* At tolerance zero, only skip exact dupes */
+					if (memcmp((char*)pt, (char*)last, pt_size) == 0)
+						continue;
+				}
 			}
-			else
-			{
-				/* exact duplicate? */
-				p1 = getPoint_internal(in, ipn-1);
-				p2 = getPoint_internal(in, ipn);
-				if (memcmp(p1, p2, ptsize) == 0)
-					/* then skip it */
-					continue;
-			}
 		}
-		/*
-		* The point is different (see above) from
-		* the previous point, so add it to output
-		*/
-		p1 = getPoint_internal(out, opn++);
-		p2 = getPoint_internal(in, ipn);
-		memcpy(p1, p2, ptsize);
 
-		prev_point = this_point;
-	}
+		/* Got to last point, and it's not very different from */
+		/* the point that preceded it. We want to keep the last */
+		/* point, not the second-to-last one, so we pull our write */
+		/* index back one value */
+		if (last_point && tolerance > 0.0 && dsq <= tolsq)
+		{
+			n_points_out--;
+		}
 
-	/* Keep the last point */
-	p1 = getPoint_internal(out, opn-1); /* Last output point */
-	p2 = getPoint_internal(in, ipn-1); /* Last input point */
-	if ( memcmp(p1, p2, ptsize) != 0 )
-		memcpy(p1, p2, ptsize);
-
-	LWDEBUGF(3, " in:%d out:%d", out->npoints, opn);
-	out->npoints = opn;
-
-	return out;
+		/* Compact all remaining values to front of array */
+		ptarray_copy_point(pa, i, n_points_out++);
+		last = pt;
+	}
+	/* Adjust array length */
+	pa->npoints = n_points_out;
+	return;
 }
 
-POINTARRAY *
-ptarray_remove_repeated_points(const POINTARRAY *in, double tolerance)
-{
-	return ptarray_remove_repeated_points_minpoints(in, tolerance, 2);
-}
 
 /************************************************************************/
 



More information about the postgis-tickets mailing list