From: Paul Ramsey <pramsey@cleverelephant.ca>
Date: Fri, 6 Mar 2015 22:03:33 +0000 (+0000)
Subject: #3074, add in lwgeom backend for subdivision
X-Git-Tag: 2.2.0rc1~607
X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=7293d5e2b5246df3a5b5c180366873b31b836bdf;p=postgis

#3074, add in lwgeom backend for subdivision


git-svn-id: http://svn.osgeo.org/postgis/trunk@13323 b70326c6-7e19-0410-871a-916f4a2858ee
---

diff --git a/liblwgeom/cunit/cu_geos.c b/liblwgeom/cunit/cu_geos.c
index 877b7bb51..8986ce931 100644
--- a/liblwgeom/cunit/cu_geos.c
+++ b/liblwgeom/cunit/cu_geos.c
@@ -18,6 +18,8 @@
 #include "lwgeom_geos.h"
 #include "cu_tester.h"
 
+#include "liblwgeom_internal.h"
+
 static void test_geos_noop(void)
 {
 	int i;
@@ -65,6 +67,38 @@ static void test_geos_noop(void)
 }
 
 
+static void test_geos_subdivide(void)
+{
+#if POSTGIS_GEOS_VERSION < 35
+	printf("%d\n", POSTGIS_GEOS_VERSION);
+	return;
+#else
+	char *ewkt = "MULTILINESTRING((0 0, 0 100))";
+	char *out_ewkt;
+	LWGEOM *geom1 = lwgeom_from_wkt(ewkt, LW_PARSER_CHECK_NONE);
+
+	LWGEOM *geom2 = lwgeom_segmentize2d(geom1, 1.0);
+	LWCOLLECTION *geom3 = lwgeom_subdivide(geom2, 80);
+	out_ewkt = lwgeom_to_ewkt((LWGEOM*)geom3);
+	// printf("\n--------\n%s\n--------\n", out_ewkt);
+	CU_ASSERT_EQUAL(2, geom3->ngeoms);
+	lwfree(out_ewkt);
+	lwcollection_free(geom3);
+	lwgeom_free(geom2);
+
+	geom2 = lwgeom_segmentize2d(geom1, 1.0);
+	geom3 = lwgeom_subdivide(geom2, 20);
+	out_ewkt = lwgeom_to_ewkt((LWGEOM*)geom3);
+	// printf("\n--------\n%s\n--------\n", out_ewkt);
+	CU_ASSERT_EQUAL(8, geom3->ngeoms);
+	lwfree(out_ewkt);
+	lwcollection_free(geom3);
+	lwgeom_free(geom2);
+
+	lwgeom_free(geom1);
+#endif
+}
+
 /*
 ** Used by test harness to register the tests in this file.
 */
@@ -73,4 +107,5 @@ void geos_suite_setup(void)
 {
 	CU_pSuite suite = CU_add_suite("GEOS", NULL, NULL);
 	PG_ADD_TEST(suite, test_geos_noop);
+	PG_ADD_TEST(suite, test_geos_subdivide);
 }
diff --git a/liblwgeom/liblwgeom.h.in b/liblwgeom/liblwgeom.h.in
index cf7c52afe..44fe32b47 100644
--- a/liblwgeom/liblwgeom.h.in
+++ b/liblwgeom/liblwgeom.h.in
@@ -2014,6 +2014,7 @@ LWGEOM *lwgeom_difference(const LWGEOM *geom1, const LWGEOM *geom2);
 LWGEOM *lwgeom_symdifference(const LWGEOM* geom1, const LWGEOM* geom2);
 LWGEOM *lwgeom_union(const LWGEOM *geom1, const LWGEOM *geom2);
 LWGEOM *lwgeom_clip_by_rect(const LWGEOM *geom1, double x0, double y0, double x1, double y1);
+LWCOLLECTION *lwgeom_subdivide(const LWGEOM *geom, int maxvertices);
 
 /**
  * Snap vertices and segments of a geometry to another using a given tolerance.
diff --git a/liblwgeom/lwgeom.c b/liblwgeom/lwgeom.c
index 186e8283c..e53c96acf 100644
--- a/liblwgeom/lwgeom.c
+++ b/liblwgeom/lwgeom.c
@@ -1872,3 +1872,103 @@ lwgeom_npoints_in_rect(const LWGEOM *geom, const GBOX *gbox)
 	return 0;
 }
 
+/* Prototype for recursion */
+static int lwgeom_subdivide_recursive(const LWGEOM *geom, int maxvertices, LWCOLLECTION *col, GBOX *clip);
+
+static int
+lwgeom_subdivide_recursive(const LWGEOM *geom, int maxvertices, LWCOLLECTION *col, GBOX *clip)
+{
+
+	/* Always just recurse into collections */
+	if ( lwgeom_is_collection(geom) )
+	{
+		LWCOLLECTION *gcol = (LWCOLLECTION*)geom;
+		int i, n = 0;
+		for ( i = 0; i < gcol->ngeoms; i++ )
+		{
+			n += lwgeom_subdivide_recursive(gcol->geoms[i], maxvertices, col, clip);
+		}
+		return n;
+	}
+	
+	/* Points and single-point multipoints can be added right into the output result */
+	if ( (geom->type == POINTTYPE) || 
+	     (geom->type == MULTIPOINTTYPE && lwgeom_count_vertices(geom) == 1) )
+	{
+		lwcollection_add_lwgeom(col, lwgeom_clone_deep(geom));
+		return 1;
+	}
+	
+	/* Haven't set up a clipping box yet? Make one then recurse again */
+	if ( ! clip )
+	{
+		double height, width;
+		const GBOX *box = lwgeom_get_bbox(geom);
+		GBOX square;
+		if ( ! box ) return 0;
+		width = box->xmax - box->xmin;
+		height = box->ymax - box->ymin;
+		if ( width > height )
+		{
+			square.xmin = box->xmin;
+			square.xmax = box->xmax;
+			square.ymin = (box->ymin + box->ymax - width)/2;
+			square.ymax = (box->ymin + box->ymax + width)/2;
+		}
+		else
+		{
+			square.ymin = box->ymin;
+			square.ymax = box->ymax;
+			square.xmin = (box->xmin + box->xmax - height)/2;
+			square.xmax = (box->xmin + box->xmax + height)/2;
+		}
+		return lwgeom_subdivide_recursive(geom, maxvertices, col, &square);
+	}
+	/* Everything is in place! Let's start subdividing! */
+	else
+	{
+		int npoints = lwgeom_npoints_in_rect(geom, clip);
+		/* Clipping rectangle didn't hit anything! No-op! */
+		if ( npoints == 0 )
+		{
+			return 0;
+		}
+		/* Clipping rectangle has reached target size: apply it! */
+		else if ( npoints < maxvertices )
+		{
+			lwcollection_add_lwgeom(col, lwgeom_clip_by_rect(geom, clip->xmin, clip->ymin, clip->xmax, clip->ymax));
+			return npoints;
+		}
+		/* Clipping rectangle too small: subdivide more! */
+		else
+		{
+			int i, n = 0;
+			GBOX boxes[4];
+			double width = FP_TOLERANCE + (clip->xmax - clip->xmin)/2;
+			for( i = 0; i < 4; i++ )
+			{
+				int ix = i / 2;
+				int iy = i % 2;
+				boxes[i].xmin = clip->xmin + ix * width;
+				boxes[i].xmax = clip->xmin + width + ix * width;
+				boxes[i].ymin = clip->ymin + iy * width;
+				boxes[i].ymax = clip->ymin + width + iy * width;
+				n += lwgeom_subdivide_recursive(geom, maxvertices, col, &(boxes[i]));
+			}
+			return n;
+		}
+	}
+	return 0;
+}
+
+LWCOLLECTION *
+lwgeom_subdivide(const LWGEOM *geom, int maxvertices)
+{
+	int n = 0;
+	LWCOLLECTION *col;
+	col = lwcollection_construct_empty(COLLECTIONTYPE, geom->srid, lwgeom_has_z(geom), lwgeom_has_m(geom));
+	n = lwgeom_subdivide_recursive(geom, maxvertices, col, NULL);
+	return col;
+}
+
+