From: Paul Ramsey Date: Sat, 7 Mar 2015 01:35:12 +0000 (+0000) Subject: #3074, subdivide any geometry into a simple set of X-Git-Tag: 2.2.0rc1~606 X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=2c31d1fc616bf4df0bfd6e64c3ae6a36925f2332;p=postgis #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 9495d4a34..fff65278c 100644 --- 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 e53c96acf..fa17c31ac 100644 --- 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 d11aa3349..de219908c 100644 --- 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 d322ecb57..44511bd8a 100644 --- 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 5518c60ca..bd52d2a6d 100644 --- 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 --------------------------------------------------------------------------------