#3074, subdivide any geometry into a simple set of

smaller parts, great for subdividing tables of
Really Big Things into things that are smaller than
the page size

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

diff --git a/liblwgeom/liblwgeom_internal.h b/liblwgeom/liblwgeom_internal.h
index 9495d4a34f04071f05ecd64cb5046bac9c125b26..fff65278c362363dafe144c2abfd2e9dcffd443f 100644 (file)
--- a/liblwgeom/liblwgeom_internal.h
+++ b/liblwgeom/liblwgeom_internal.h
@@ -439,3 +439,4 @@ extern int _lwgeom_interrupt_requested;
 int ptarray_npoints_in_rect(const POINTARRAY *pa, const GBOX *gbox);
 int gbox_contains_point2d(const GBOX *g, const POINT2D *p);
 int lwgeom_npoints_in_rect(const LWGEOM *geom, const GBOX *gbox);
+int lwpoly_contains_point(const LWPOLY *poly, const POINT2D *pt);

diff --git a/liblwgeom/lwgeom.c b/liblwgeom/lwgeom.c
index e53c96acf676f4cee9d729d99e1cd0588862da53..fa17c31acf1f33cf8fbb6f9f24cfc0d9e07b76de 100644 (file)
--- a/liblwgeom/lwgeom.c
+++ b/liblwgeom/lwgeom.c
@@ -1905,6 +1905,7 @@ lwgeom_subdivide_recursive(const LWGEOM *geom, int maxvertices, LWCOLLECTION *co
                double height, width;
                const GBOX *box = lwgeom_get_bbox(geom);
                GBOX square;
+               gbox_init(&square);
                if ( ! box ) return 0;
                width = box->xmax - box->xmin;
                height = box->ymax - box->ymin;
@@ -1928,16 +1929,65 @@ lwgeom_subdivide_recursive(const LWGEOM *geom, int maxvertices, LWCOLLECTION *co
        else
        {
                int npoints = lwgeom_npoints_in_rect(geom, clip);
-               /* Clipping rectangle didn't hit anything! No-op! */
+               /* Clipping rectangle didn't hit anything! */
                if ( npoints == 0 )
                {
-                       return 0;
+                       /* Is it because it's a rectangle totally *inside* a polygon? */
+                       if ( geom->type == POLYGONTYPE )
+                       {
+                               POINT2D pt;
+                               pt.x = clip->xmin;
+                               pt.y = clip->ymin;
+                               if ( lwpoly_contains_point((LWPOLY*)geom, &pt) )
+                               {
+                                       LWGEOM *clipped = lwgeom_clip_by_rect(geom, clip->xmin, clip->ymin, clip->xmax, clip->ymax);
+                                       /* Hm, clipping left nothing behind, skip it */
+                                       if ( lwgeom_is_empty(clipped) )
+                                       {
+                                               return 0;
+                                       }
+                                       /* Add the clipped part */
+                                       else
+                                       {
+                                               lwcollection_add_lwgeom(col, clipped);
+                                               return 5;
+                                       }
+                               }
+                               else
+                               {
+                                       return 0;
+                               }
+                       }
+                       else
+                       {
+                               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;
+                       /* Oh, it also covers the whole geometry */
+                       if ( npoints == lwgeom_count_vertices(geom) )
+                       {
+                               lwcollection_add_lwgeom(col, lwgeom_clone_deep(geom));
+                               return npoints;
+                       }
+                       /* Partial coverage, we need to clip */
+                       else
+                       {
+                               LWGEOM *clipped = lwgeom_clip_by_rect(geom, clip->xmin, clip->ymin, clip->xmax, clip->ymax);
+                               /* Hm, clipping left nothing behind, skip it */
+                               if ( lwgeom_is_empty(clipped) )
+                               {
+                                       return 0;
+                               }
+                               /* Add the clipped part */
+                               else
+                               {
+                                       lwcollection_add_lwgeom(col, clipped);
+                                       return npoints;
+                               }                               
+                       }                               
                }
                /* Clipping rectangle too small: subdivide more! */
                else
@@ -1949,6 +1999,7 @@ lwgeom_subdivide_recursive(const LWGEOM *geom, int maxvertices, LWCOLLECTION *co
                        {
                                int ix = i / 2;
                                int iy = i % 2;
+                               gbox_init(&(boxes[i]));
                                boxes[i].xmin = clip->xmin + ix * width;
                                boxes[i].xmax = clip->xmin + width + ix * width;
                                boxes[i].ymin = clip->ymin + iy * width;

diff --git a/liblwgeom/lwpoly.c b/liblwgeom/lwpoly.c
index d11aa3349822669568fa2f859868605ae8373c2f..de219908cc95164b8232ee4e4abd369162f303db 100644 (file)
--- a/liblwgeom/lwpoly.c
+++ b/liblwgeom/lwpoly.c
@@ -489,6 +489,26 @@ lwpoly_startpoint(const LWPOLY* poly, POINT4D* pt)
        return ptarray_startpoint(poly->rings[0], pt);
 }
 
+int
+lwpoly_contains_point(const LWPOLY *poly, const POINT2D *pt)
+{
+       int i;
+       
+       if ( lwpoly_is_empty(poly) )
+               return LW_FALSE;
+       
+       if ( ptarray_contains_point(poly->rings[0], pt) == LW_OUTSIDE )
+               return LW_FALSE;
+       
+       for ( i = 1; i < poly->nrings; i++ )
+       {
+               if ( ptarray_contains_point(poly->rings[i], pt) == LW_INSIDE )
+                       return LW_FALSE;
+       }
+       return LW_TRUE;
+}
+
+
 
 LWPOLY* lwpoly_grid(const LWPOLY *poly, const gridspec *grid)
 {

diff --git a/postgis/lwgeom_dump.c b/postgis/lwgeom_dump.c
index d322ecb57fdbf85ac6cabcf411e49a708b910973..44511bd8afb9548fe080a2c44efff11d303e96d7 100644 (file)
--- a/postgis/lwgeom_dump.c
+++ b/postgis/lwgeom_dump.c
@@ -31,6 +31,7 @@
 
 Datum LWGEOM_dump(PG_FUNCTION_ARGS);
 Datum LWGEOM_dump_rings(PG_FUNCTION_ARGS);
+Datum ST_Subdivide(PG_FUNCTION_ARGS);
 
 typedef struct GEOMDUMPNODE_T
 {
@@ -300,3 +301,104 @@ Datum LWGEOM_dump_rings(PG_FUNCTION_ARGS)
 
 }
 
+
+struct FLATCOLLECTIONDUMPSTATE
+{
+       int geomnum;
+       LWCOLLECTION *col;
+};
+
+/*
+* Break an object up into smaller objects of no more than N vertices 
+*/
+PG_FUNCTION_INFO_V1(ST_Subdivide);
+Datum ST_Subdivide(PG_FUNCTION_ARGS)
+{
+#if POSTGIS_GEOS_VERSION < 35
+
+       elog(ERROR, ("The GEOS version this PostGIS binary "
+               "was compiled against (%d) doesn't support "
+               "'%s' function (3.5.0+ required)",
+               POSTGIS_GEOS_VERSION, __func__);
+       PG_RETURN_NULL();
+
+#else /* POSTGIS_GEOS_VERSION >= 35 */
+
+       typedef struct
+       {
+               int nextgeom;
+               int numgeoms;
+               LWCOLLECTION *col;
+       } collection_fctx;
+
+       FuncCallContext *funcctx;
+       collection_fctx *fctx;
+       MemoryContext oldcontext;
+
+       /* stuff done only on the first call of the function */
+       if (SRF_IS_FIRSTCALL())
+       {
+               GSERIALIZED *gser;
+               LWGEOM *geom;
+               LWCOLLECTION *col;
+               int maxvertices = 256;
+
+               /* create a function context for cross-call persistence */
+               funcctx = SRF_FIRSTCALL_INIT();
+
+               /*
+               * switch to memory context appropriate for multiple function calls
+               */
+               oldcontext = MemoryContextSwitchTo(funcctx->multi_call_memory_ctx);
+
+               /*
+               * Get the geometry value 
+               */
+               gser = PG_GETARG_GSERIALIZED_P(0);
+               geom = lwgeom_from_gserialized(gser);
+               
+               /*
+               * Get the max vertices value 
+               */
+               maxvertices = PG_GETARG_INT32(1);
+               
+               /*
+               * Compute the subdivision of the geometry
+               */
+               col = lwgeom_subdivide(geom, maxvertices);
+               
+               if ( ! col ) 
+                       SRF_RETURN_DONE(funcctx);
+
+               /* allocate memory for user context */
+               fctx = (collection_fctx *) palloc(sizeof(collection_fctx));
+
+               /* initialize state */
+               fctx->nextgeom = 0;
+               fctx->numgeoms = col->ngeoms;
+               fctx->col = col;
+
+               /* save user context, switch back to function context */
+               funcctx->user_fctx = fctx;
+               MemoryContextSwitchTo(oldcontext);
+       }
+
+       /* stuff done on every call of the function */
+       funcctx = SRF_PERCALL_SETUP();
+       fctx = funcctx->user_fctx;
+
+       if (fctx->nextgeom < fctx->numgeoms)
+       {
+               GSERIALIZED *gpart = geometry_serialize(fctx->col->geoms[fctx->nextgeom]);
+               fctx->nextgeom++;
+               SRF_RETURN_NEXT(funcctx, PointerGetDatum(gpart));
+       }
+       else
+       {
+               /* do when there is no more left */
+               SRF_RETURN_DONE(funcctx);
+       }
+       
+#endif /* POSTGIS_GEOS_VERSION >= 35 */
+}
+

diff --git a/postgis/postgis.sql.in b/postgis/postgis.sql.in
index 5518c60caa268b67d878c8dd8713faac3fcddd9e..bd52d2a6d708a4f0190a0c3041665f9ab9e6cae3 100644 (file)
--- a/postgis/postgis.sql.in
+++ b/postgis/postgis.sql.in
@@ -3221,6 +3221,15 @@ CREATE OR REPLACE FUNCTION ST_ClipByBox2d(geom geometry, box box2d)
        LANGUAGE 'c' IMMUTABLE STRICT
        COST 50;
 
+-- Requires GEOS >= 3.5.0
+-- Availability: 2.2.0
+CREATE OR REPLACE FUNCTION ST_Subdivide(geom geometry, maxvertices integer DEFAULT 256)
+       RETURNS setof geometry
+       AS 'MODULE_PATHNAME', 'ST_Subdivide'
+       LANGUAGE 'c' IMMUTABLE STRICT
+       COST 100;
+
+
 --------------------------------------------------------------------------------
 -- ST_CleanGeometry / ST_MakeValid
 --------------------------------------------------------------------------------