ST_Subdivide: improved handling of polygons with many holes.

Closes #4038
Closes #202



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

NEWS

@@ -46,6 +46,9 @@ PostGIS 2.5.0
   - #4037, Invalid input geometry is fixed with MakeValid for GEOS exceptions in
            ST_Intersection, ST_Union, ST_Difference, ST_SymDifference (Darafei
            Praliaskouski)
+  - #4038, ST_Subdivide now selects pivot for geometry split that reuses input
+           vertices. ST_ClipByBox2D is stubbed with ST_Intersection because of
+           robustness issues. (Darafei Praliaskouski)
 
 PostGIS 2.4.0
 2017/09/30

liblwgeom/cunit/cu_geos.c

@@ -112,13 +112,13 @@ test_geos_offsetcurve(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 *ewkt = "LINESTRING(0 0, 10 10)";
 	char *out_ewkt;
 	LWGEOM *geom1 = lwgeom_from_wkt(ewkt, LW_PARSER_CHECK_NONE);
-	LWGEOM *geom2 = lwgeom_segmentize2d(geom1, 1.0);
+	/* segmentize as geography to generate a non-simple curve */
+	LWGEOM *geom2 = lwgeom_segmentize_sphere(geom1, 0.002);
 
 	LWCOLLECTION *geom3 = lwgeom_subdivide(geom2, 80);
 	out_ewkt = lwgeom_to_ewkt((LWGEOM*)geom3);

liblwgeom/lwgeom.c

@@ -2248,24 +2248,33 @@ lwgeom_grid(const LWGEOM *lwgeom, const gridspec *grid)
 
 
 /* Prototype for recursion */
-static int
-lwgeom_subdivide_recursive(const LWGEOM *geom, uint32_t maxvertices, uint32_t depth, LWCOLLECTION *col, const GBOX *clip);
+static int lwgeom_subdivide_recursive(const LWGEOM *geom, uint32_t maxvertices, uint32_t depth, LWCOLLECTION *col);
 
 static int
-lwgeom_subdivide_recursive(const LWGEOM *geom, uint32_t maxvertices, uint32_t depth, LWCOLLECTION *col, const GBOX *clip)
+lwgeom_subdivide_recursive(const LWGEOM *geom, uint32_t maxvertices, uint32_t depth, LWCOLLECTION *col)
 {
 	const uint32_t maxdepth = 50;
+	GBOX *clip = gbox_copy(lwgeom_get_bbox(geom));
 	uint32_t nvertices = 0;
 	uint32_t i, n = 0;
-	double width = clip->xmax - clip->xmin;
-	double height = clip->ymax - clip->ymin;
-	GBOX subbox1, subbox2;
-	LWGEOM *clipped1, *clipped2;
+	uint32_t split_ordinate;
+	double width;
+	double height;
+	double pivot = DBL_MAX;
+	double center = DBL_MAX;
+	LWPOLY *lwpoly = NULL;
+
+	GBOX *subbox1;
+	GBOX *subbox2;
+	LWGEOM *clipped;
+
+	if (!clip) return 0;
+
+	width = clip->xmax - clip->xmin;
+	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 )
 	{
@@ -2275,21 +2284,31 @@ lwgeom_subdivide_recursive(const LWGEOM *geom, uint32_t maxvertices, uint32_t de
 			return 1;
 		}
 		else
-		{
 			return 0;
-		}
+	}
+
+	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;
 		int n = 0;
+		/* Don't increment depth yet, since we aren't actually
+		 * subdividing geomtries yet */
 		for ( i = 0; i < incol->ngeoms; i++ )
-		{
-			/* Don't increment depth yet, since we aren't actually subdividing geomtries yet */
-			n += lwgeom_subdivide_recursive(incol->geoms[i], maxvertices, depth, col, clip);
-		}
+			n += lwgeom_subdivide_recursive(incol->geoms[i], maxvertices, depth, col);
 		return n;
 	}
 
@@ -2302,73 +2321,102 @@ lwgeom_subdivide_recursive(const LWGEOM *geom, uint32_t maxvertices, uint32_t de
 	}
 
 	nvertices = lwgeom_count_vertices(geom);
+
 	/* Skip empties entirely */
-	if ( nvertices == 0 )
-	{
-		return 0;
-	}
+	if (nvertices == 0) return 0;
 
 	/* If it is under the vertex tolerance, just add it, we're done */
-	if ( nvertices < maxvertices )
+	if (nvertices <= maxvertices)
 	{
 		lwcollection_add_lwgeom(col, lwgeom_clone_deep(geom));
 		return 1;
 	}
 
-	subbox1 = subbox2 = *clip;
-	if ( width > height )
-	{
-		subbox1.xmax = subbox2.xmin = (clip->xmin + clip->xmax)/2;
-	}
+	split_ordinate = (width > height) ? 0 : 1;
+	if (split_ordinate == 0)
+		center = (clip->xmin + clip->xmax) / 2;
 	else
-	{
-		subbox1.ymax = subbox2.ymin = (clip->ymin + clip->ymax)/2;
-	}
+		center = (clip->ymin + clip->ymax) / 2;
 
-	if ( height == 0 )
+	if (geom->type == POLYGONTYPE)
 	{
-		subbox1.ymax += FP_TOLERANCE;
-		subbox2.ymax += FP_TOLERANCE;
-		subbox1.ymin -= FP_TOLERANCE;
-		subbox2.ymin -= FP_TOLERANCE;
-	}
+		uint32_t ring_to_trim = 0;
+		double ring_area = 0;
+		double pivot_eps = DBL_MAX;
+		double pt_eps = DBL_MAX;
+		POINTARRAY *pa;
+		pivot = DBL_MAX;
+		lwpoly = (LWPOLY *)geom;
 
-	if ( width == 0 )
-	{
-		subbox1.xmax += FP_TOLERANCE;
-		subbox2.xmax += FP_TOLERANCE;
-		subbox1.xmin -= FP_TOLERANCE;
-		subbox2.xmin -= FP_TOLERANCE;
+		/* if there are more points in holes than in outer ring */
+		if (nvertices > 2 * lwpoly->rings[0]->npoints)
+		{
+			/* trim holes starting from biggest */
+			for (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 (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;
+			}
+		}
 	}
 
-	clipped1 = lwgeom_clip_by_rect(geom, subbox1.xmin, subbox1.ymin, subbox1.xmax, subbox1.ymax);
-	clipped2 = lwgeom_clip_by_rect(geom, subbox2.xmin, subbox2.ymin, subbox2.xmax, subbox2.ymax);
+	subbox1 = gbox_copy(clip);
+	subbox2 = gbox_copy(clip);
+
+	if (pivot == DBL_MAX) pivot = center;
+
+	if (split_ordinate == 0)
+		subbox1->xmax = subbox2->xmin = pivot;
+	else
+		subbox1->ymax = subbox2->ymin = pivot;
 
 	++depth;
 
-	if ( clipped1 )
+	clipped = lwgeom_clip_by_rect(geom, subbox1->xmin, subbox1->ymin, subbox1->xmax, subbox1->ymax);
+	if (clipped)
 	{
-		n += lwgeom_subdivide_recursive(clipped1, maxvertices, depth, col, &subbox1);
-		lwgeom_free(clipped1);
+		n += lwgeom_subdivide_recursive(clipped, maxvertices, depth, col);
+		lwgeom_free(clipped);
 	}
 
-	if ( clipped2 )
+	clipped = lwgeom_clip_by_rect(geom, subbox2->xmin, subbox2->ymin, subbox2->xmax, subbox2->ymax);
+	if (clipped)
 	{
-		n += lwgeom_subdivide_recursive(clipped2, maxvertices, depth, col, &subbox2);
-		lwgeom_free(clipped2);
+		n += lwgeom_subdivide_recursive(clipped, maxvertices, depth, col);
+		lwgeom_free(clipped);
 	}
 
 	return n;
-
 }
 
 LWCOLLECTION *
 lwgeom_subdivide(const LWGEOM *geom, uint32_t maxvertices)
 {
 	static uint32_t startdepth = 0;
-	static uint32_t minmaxvertices = 8;
+	static uint32_t minmaxvertices = 5;
 	LWCOLLECTION *col;
-	GBOX clip;
 
 	col = lwcollection_construct_empty(COLLECTIONTYPE, geom->srid, lwgeom_has_z(geom), lwgeom_has_m(geom));
 
@@ -2381,8 +2429,7 @@ lwgeom_subdivide(const LWGEOM *geom, uint32_t maxvertices)
 		lwerror("%s: cannot subdivide to fewer than %d vertices per output", __func__, minmaxvertices);
 	}
 
-	clip = *(lwgeom_get_bbox(geom));
-	lwgeom_subdivide_recursive(geom, maxvertices, startdepth, col, &clip);
+	lwgeom_subdivide_recursive(geom, maxvertices, startdepth, col);
 	lwgeom_set_srid((LWGEOM*)col, geom->srid);
 	return col;
 }

liblwgeom/lwgeom_geos.c

@@ -315,6 +315,13 @@ LWGEOM2GEOS(const LWGEOM* lwgeom, uint8_t autofix)
 	char* wkt;
 #endif
 
+	if (autofix)
+	{
+		/* cross fingers and try without autofix, maybe it'll work? */
+		g = LWGEOM2GEOS(lwgeom, LW_FALSE);
+		if (g) return g;
+	}
+
 	LWDEBUGF(4, "LWGEOM2GEOS got a %s", lwtype_name(lwgeom->type));
 
 	if (lwgeom_has_arc(lwgeom))
@@ -872,17 +879,35 @@ lwgeom_union(const LWGEOM* geom1, const LWGEOM* geom2)
 	return result;
 }
 
-LWGEOM*
-lwgeom_clip_by_rect(const LWGEOM* geom, double x0, double y0, double x1, double y1)
+LWGEOM *
+lwgeom_clip_by_rect(const LWGEOM *geom1, double x1, double y1, double x2, double y2)
 {
-#if POSTGIS_GEOS_VERSION < 35
-	lwerror(
-	    "The GEOS version this postgis binary was compiled against (%d) doesn't support 'GEOSClipByRect' function "
-	    "(3.5.0+ required)",
-	    POSTGIS_GEOS_VERSION);
-	return NULL;
-#else  /* POSTGIS_GEOS_VERSION >= 35 */
 	LWGEOM* result;
+	LWGEOM *tmp;
+
+	result = lwgeom_intersection(geom1, (LWGEOM *)lwpoly_construct_envelope(geom1->srid, x1, y1, x2, y2));
+
+	if (!result) return NULL;
+
+	/* clipping should not produce lower dimension objects */
+	if (
+	    /* input has exact dimensionality, isn't a generic collection */
+	    geom1->type != COLLECTIONTYPE &&
+	    /* output may have different things inside */
+	    result->type == COLLECTIONTYPE)
+	{
+		tmp = lwcollection_as_lwgeom(lwcollection_extract(lwgeom_as_lwcollection(result), lwgeom_dimension(geom1) + 1));
+		lwfree(result);
+		result = tmp;
+		if (!result) return NULL;
+	}
+
+	/* clean up stray points on geometry boundary */
+	lwgeom_simplify_in_place(result, 0.0, LW_TRUE);
+
+	return result;
+
+#if 0  /* POSTGIS_GEOS_VERSION >= 35, enable only after bugs in geos are fixed */
 	int32_t srid = get_result_srid(geom, NULL, __func__);
 	uint8_t is3d = FLAGS_GET_Z(geom->flags);
 	GEOSGeometry *g1, *g3;
@@ -896,7 +921,7 @@ lwgeom_clip_by_rect(const LWGEOM* geom, double x0, double y0, double x1, double
 
 	if (!input_lwgeom_to_geos(&g1, geom, __func__)) return NULL;
 
-	g3 = GEOSClipByRect(g1, x0, y0, x1, y1);
+	g3 = GEOSClipByRect(g1, x1, y1, x2, y2);
 
 	if (!g3) return geos_clean_and_fail(g1, NULL, NULL, __func__);
 

regress/clipbybox2d.sql

@@ -1,28 +1,24 @@
 -- Overlap
-SELECT ST_ClipByBox2d(ST_MakeEnvelope(0,0,10,10),ST_MakeEnvelope(5,5,20,20))::box2d;
+SELECT '1', ST_ClipByBox2d(ST_MakeEnvelope(0,0,10,10),ST_MakeEnvelope(5,5,20,20))::box2d;
 -- Geom covers rect
-SELECT ST_ClipByBox2d(ST_MakeEnvelope(0,0,10,10),ST_MakeEnvelope(5,5,8,8))::box2d;
+SELECT '2', ST_ClipByBox2d(ST_MakeEnvelope(0,0,10,10),ST_MakeEnvelope(5,5,8,8))::box2d;
 -- Geom within rect
-SELECT ST_ClipByBox2d(ST_MakeEnvelope(2,2,8,8),ST_MakeEnvelope(0,0,10,10))::box2d;
+SELECT '3', ST_ClipByBox2d(ST_MakeEnvelope(2,2,8,8),ST_MakeEnvelope(0,0,10,10))::box2d;
 -- Multipoint with point inside, point outside and point on boundary
-SELECT ST_AsText(ST_ClipByBox2d('MULTIPOINT(-1 -1, 0 0, 2 2)'::geometry,ST_MakeEnvelope(0,0,10,10)));
--- Invalid polygon (bow-tie) -- ST_Intersection throws an exception here
-SELECT ST_AsText(ST_ClipByBox2d('POLYGON((0 0, 10 10, 0 10, 10 0, 0 0))', ST_MakeEnvelope(2,2,8,8)));
+SELECT '4', ST_AsText(ST_ClipByBox2d('MULTIPOINT(-1 -1, 0 0, 2 2)'::geometry,ST_MakeEnvelope(0,0,10,10)));
+-- Invalid polygon (bow-tie)
+SELECT '5', ST_AsText(ST_ClipByBox2d('POLYGON((0 0, 10 10, 0 10, 10 0, 0 0))', ST_MakeEnvelope(2,2,8,8)));
 -- Invalid polygon (lineal self-intersection) -- ST_Intersection returns a collection
-SELECT ST_AsText(ST_ClipByBox2d('POLYGON((0 0,5 4,5 6,0 10,10 10,5 6,5 4,10 0,0 0))', ST_MakeEnvelope(2,2,10,5)));
+SELECT '6', ST_AsText(ST_ClipByBox2d('POLYGON((0 0,5 4,5 6,0 10,10 10,5 6,5 4,10 0,0 0))', ST_MakeEnvelope(2,2,10,5)));
 -- Invalid polygon (non-closed ring) -- ST_Intersection raises an exception
 -- The HEXWKB represents 'POLYGON((0 0, 10 0, 10 10, 0 10))'
-SELECT ST_AsText(ST_ClipByBox2d(
+SELECT '7', ST_AsText(ST_ClipByBox2d(
   '0103000000010000000400000000000000000000000000000000000000000000000000244000000000000000000000000000002440000000000000244000000000000000000000000000002440'
   , ST_MakeEnvelope(2,2,5,5)));
 -- Geometry disjoint from box, from a table
 -- See http://trac.osgeo.org/postgis/ticket/2950
 CREATE TEMPORARY TABLE t AS SELECT
-'SRID=3857;POLYGON((41 20,41 0,21 0,1 20,1 40,21 40,41 20))'
-::geometry g;
-SELECT ST_AsEWKT(ST_ClipByBox2d(g, ST_MakeEnvelope(-20,-20,-10,-10))) FROM t;
+'SRID=3857;POLYGON((41 20,41 0,21 0,1 20,1 40,21 40,41 20))'::geometry g;
+SELECT '8', ST_AsEWKT(ST_ClipByBox2d(g, ST_MakeEnvelope(-20,-20,-10,-10))) FROM t;
 -- See http://trac.osgeo.org/postgis/ticket/2954
-SELECT ST_AsEWKT(ST_ClipByBox2D('SRID=4326;POINT(0 0)','BOX3D(-1 -1,1 1)'::box3d::box2d));
-
-SELECT '#3135', st_astext(ST_SubDivide(ST_GeomFromText('POLYGON((1 2,1 2,1 2,1 2))'), 2));
-
+SELECT '9', ST_AsEWKT(ST_ClipByBox2D('SRID=4326;POINT(0 0)','BOX3D(-1 -1,1 1)'::box3d::box2d));

regress/clipbybox2d_expected

@@ -1,10 +1,13 @@
-BOX(5 5,10 10)
-BOX(5 5,8 8)
-BOX(2 2,8 8)
-POINT(2 2)
-POLYGON((2 2,8 2,2 8,8 8,2 2))
-POLYGON((2.5 2,5 4,5 5,5 4,7.5 2,2.5 2))
-POLYGON((2 2,2 5,5 5,5 2,2 2))
-SRID=3857;POLYGON EMPTY
-SRID=4326;POINT(0 0)
-ERROR:  lwgeom_subdivide: cannot subdivide to fewer than 8 vertices per output
+1|BOX(5 5,10 10)
+2|BOX(5 5,8 8)
+3|BOX(2 2,8 8)
+4|MULTIPOINT(0 0,2 2)
+NOTICE:  lwgeom_intersection: GEOS Error: TopologyException: Input geom 0 is invalid: Self-intersection
+NOTICE:  Self-intersection
+NOTICE:  Your geometry dataset is not valid per OGC Specification. Please fix it with manual review of entries that are not ST_IsValid(geom). Retrying GEOS operation with ST_MakeValid of your input.
+NOTICE:  Self-intersection
+5|MULTIPOLYGON(((2 2,5 5,8 2,2 2)),((5 5,2 8,8 8,5 5)))
+6|MULTIPOLYGON(((2.5 2,5 4,5 5,10 5,10 2,2.5 2)))
+7|POLYGON((2 2,2 5,5 5,5 2,2 2))
+8|SRID=3857;POLYGON EMPTY
+9|SRID=4326;POINT(0 0)

regress/regress_big_polygon.sql

@@ -0,0 +1,9 @@
+\set QUIET on
+SET client_min_messages TO WARNING;
+drop table if exists big_polygon;
+create table big_polygon(geom geometry);
+alter table big_polygon alter column geom set storage main;
+-- Data (c) OpenStreetMap contributors, ODbL / http://osm.org/
+\copy big_polygon from 'big_polygon.wkb'
+\set QUIET off
+SET client_min_messages TO NOTICE;