[postgis-tickets] r16676 - Make lwgeom_subdivide use lwgeom_intersection for clipping

[Daniel Baston]

Author: dbaston
Date: 2018-07-31 17:45:01 -0700 (Tue, 31 Jul 2018)
New Revision: 16676

Modified:
   trunk/liblwgeom/lwgeom.c
Log:
Make lwgeom_subdivide use lwgeom_intersection for clipping

lwgeom_clip_by_rect describes itself as clipping in a "fast but possibly
dirty way." This is not a good fit for lwgeom_subdivide, a primary use of
which is the optimization of point-in-polygon processes. That application
requires that the area covered by the subdivided geometry be the same as
the area covered by the original geometry, which in turn requires that a
robust intersection routine be used.

References #4038
Closes https://github.com/postgis/postgis/pull/281



Modified: trunk/liblwgeom/lwgeom.c
===================================================================
--- trunk/liblwgeom/lwgeom.c	2018-07-30 16:38:45 UTC (rev 16675)
+++ trunk/liblwgeom/lwgeom.c	2018-08-01 00:45:01 UTC (rev 16676)
@@ -2249,10 +2249,10 @@
 
 
 /* Prototype for recursion */
-static int lwgeom_subdivide_recursive(const LWGEOM *geom, uint32_t maxvertices, uint32_t depth, LWCOLLECTION *col);
+static int lwgeom_subdivide_recursive(const LWGEOM *geom, uint8_t dimension, uint32_t maxvertices, uint32_t depth, LWCOLLECTION *col);
 
 static int
-lwgeom_subdivide_recursive(const LWGEOM *geom, uint32_t maxvertices, uint32_t depth, LWCOLLECTION *col)
+lwgeom_subdivide_recursive(const LWGEOM *geom, uint8_t dimension, uint32_t maxvertices, uint32_t depth, LWCOLLECTION *col)
 {
 	const uint32_t maxdepth = 50;
 	GBOX clip, subbox1, subbox2;
@@ -2275,7 +2275,7 @@
 
 	if ( width == 0.0 && height == 0.0 )
 	{
-		if ( geom->type == POINTTYPE )
+		if ( geom->type == POINTTYPE && dimension == 0)
 		{
 			lwcollection_add_lwgeom(col, lwgeom_clone_deep(geom));
 			return 1;
@@ -2303,12 +2303,19 @@
 		LWCOLLECTION *incol = (LWCOLLECTION*)geom;
 		int n = 0;
 		/* Don't increment depth yet, since we aren't actually
-		 * subdividing geomtries yet */
+		 * subdividing geometries yet */
 		for ( i = 0; i < incol->ngeoms; i++ )
-			n += lwgeom_subdivide_recursive(incol->geoms[i], maxvertices, depth, col);
+			n += lwgeom_subdivide_recursive(incol->geoms[i], dimension, maxvertices, depth, col);
 		return n;
 	}
 
+	if (lwgeom_dimension(geom) < dimension)
+	{
+		/* We've hit a lower dimension object produced by clipping at
+		 * a shallower recursion level. Ignore it. */
+		return 0;
+	}
+
 	/* But don't go too far. 2^50 ~= 10^15, that's enough subdivision */
 	/* Just add what's left */
 	if ( depth > maxdepth )
@@ -2391,17 +2398,23 @@
 
 	++depth;
 
-	clipped = lwgeom_clip_by_rect(geom, subbox1.xmin, subbox1.ymin, subbox1.xmax, subbox1.ymax);
+	LWGEOM* subbox = (LWGEOM*) lwpoly_construct_envelope(geom->srid, subbox1.xmin, subbox1.ymin, subbox1.xmax, subbox1.ymax);
+	clipped = lwgeom_intersection(geom, subbox);
+	lwgeom_simplify_in_place(clipped, 0.0, LW_TRUE);
+	lwgeom_free(subbox);
 	if (clipped)
 	{
-		n += lwgeom_subdivide_recursive(clipped, maxvertices, depth, col);
+		n += lwgeom_subdivide_recursive(clipped, dimension, maxvertices, depth, col);
 		lwgeom_free(clipped);
 	}
 
-	clipped = lwgeom_clip_by_rect(geom, subbox2.xmin, subbox2.ymin, subbox2.xmax, subbox2.ymax);
+	subbox = (LWGEOM*) lwpoly_construct_envelope(geom->srid, subbox2.xmin, subbox2.ymin, subbox2.xmax, subbox2.ymax);
+	clipped = lwgeom_intersection(geom, subbox);
+	lwgeom_simplify_in_place(clipped, 0.0, LW_TRUE);
+	lwgeom_free(subbox);
 	if (clipped)
 	{
-		n += lwgeom_subdivide_recursive(clipped, maxvertices, depth, col);
+		n += lwgeom_subdivide_recursive(clipped, dimension, maxvertices, depth, col);
 		lwgeom_free(clipped);
 	}
 
@@ -2426,7 +2439,7 @@
 		lwerror("%s: cannot subdivide to fewer than %d vertices per output", __func__, minmaxvertices);
 	}
 
-	lwgeom_subdivide_recursive(geom, maxvertices, startdepth, col);
+	lwgeom_subdivide_recursive(geom, lwgeom_dimension(geom), maxvertices, startdepth, col);
 	lwgeom_set_srid((LWGEOM*)col, geom->srid);
 	return col;
 }