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 6038230ae079ebbfb29bf98b34cb6541b5823b3a..d582f02da180ee0ac709646ff9e8f303bbcd55c9 100644 (file)
--- 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;
 }