From: Daniel Baston <dbaston@gmail.com>
Date: Wed, 1 Aug 2018 00:45:01 +0000 (+0000)
Subject: Make lwgeom_subdivide use lwgeom_intersection for clipping
X-Git-Tag: 2.5.0beta2~14
X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=12b917fb3116ea6c328dc573ca0f5c910654bfb2;p=postgis

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



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

diff --git a/liblwgeom/lwgeom.c b/liblwgeom/lwgeom.c
index 6038230ae..d582f02da 100644
--- a/liblwgeom/lwgeom.c
+++ b/liblwgeom/lwgeom.c
@@ -2249,10 +2249,10 @@ 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);
+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 @@ lwgeom_subdivide_recursive(const LWGEOM *geom, uint32_t maxvertices, uint32_t de
 
 	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 @@ lwgeom_subdivide_recursive(const LWGEOM *geom, uint32_t maxvertices, uint32_t de
 		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 @@ lwgeom_subdivide_recursive(const LWGEOM *geom, uint32_t maxvertices, uint32_t de
 
 	++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 @@ lwgeom_subdivide(const LWGEOM *geom, uint32_t maxvertices)
 		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;
 }