From b7e1ad5dd11579cb7353bffeb67262a1801d05d5 Mon Sep 17 00:00:00 2001 From: Sandro Santilli Date: Mon, 15 Mar 2010 18:00:12 +0000 Subject: [PATCH] Implement split-poly-by-line git-svn-id: http://svn.osgeo.org/postgis/trunk@5437 b70326c6-7e19-0410-871a-916f4a2858ee --- postgis/lwgeom_geos_split.c | 163 ++++++++++++++++++++++++++++++++++-- regress/split.sql | 12 +++ regress/split_expected | 4 + 3 files changed, 173 insertions(+), 6 deletions(-) diff --git a/postgis/lwgeom_geos_split.c b/postgis/lwgeom_geos_split.c index 39698e82e..f88986f9c 100644 --- a/postgis/lwgeom_geos_split.c +++ b/postgis/lwgeom_geos_split.c @@ -56,7 +56,6 @@ lwline_split_by_line(LWLINE* lwline_in, LWLINE* blade_in) GEOSGeometry* gdiff; /* difference */ GEOSGeometry* g1; GEOSGeometry* g2; - char* i9; int ret; /* Possible outcomes: @@ -64,9 +63,7 @@ lwline_split_by_line(LWLINE* lwline_in, LWLINE* blade_in) * 1. The lines do not cross or overlap * -> Return a collection with single element * 2. The lines cross - * -> Return a collection 2 elements: - * o segments falling on the left of the directed blade - * o segments falling on the right of the directed blade + * -> Return a collection of all elements resulting from the split */ initGEOS(lwgeom_geos_error, lwgeom_geos_error); @@ -218,13 +215,167 @@ lwline_split(LWLINE* lwline_in, LWGEOM* blade_in) return NULL; } -static LWGEOM* lwpoly_split(LWPOLY* lwgeom_in, LWGEOM* blade_in); +/* Initializes and uses GEOS internally */ +static LWGEOM* lwpoly_split_by_line(LWPOLY* lwgeom_in, LWLINE* blade_in); +static LWGEOM* +lwpoly_split_by_line(LWPOLY* lwpoly_in, LWLINE* blade_in) +{ + LWCOLLECTION* out; + GEOSGeometry* g1; + GEOSGeometry* g2; + GEOSGeometry* g1_bounds; + GEOSGeometry* polygons; + const GEOSGeometry *vgeoms[1]; + int i,n; + + /* Possible outcomes: + * + * 1. The line does not split the polygon + * -> Return a collection with single element + * 2. The line does split the polygon + * -> Return a collection of all elements resulting from the split + */ + + initGEOS(lwgeom_geos_error, lwgeom_geos_error); + + g1 = LWGEOM2GEOS((LWGEOM*)lwpoly_in); + if ( NULL == g1 ) { + lwerror("LWGEOM2GEOS: %s", lwgeom_geos_errmsg); + return NULL; + } + g1_bounds = GEOSBoundary(g1); + if ( NULL == g1_bounds ) + { + GEOSGeom_destroy(g1); + lwerror("GEOSBoundary: %s", lwgeom_geos_errmsg); + return NULL; + } + + g2 = LWGEOM2GEOS((LWGEOM*)blade_in); + if ( NULL == g2 ) { + GEOSGeom_destroy(g1); + GEOSGeom_destroy(g1_bounds); + lwerror("LWGEOM2GEOS: %s", lwgeom_geos_errmsg); + return NULL; + } + + vgeoms[0] = GEOSUnion(g1_bounds, g2); + if ( NULL == vgeoms[0] ) + { + GEOSGeom_destroy(g1); + GEOSGeom_destroy(g2); + GEOSGeom_destroy(g1_bounds); + lwerror("GEOSUnion: %s", lwgeom_geos_errmsg); + return NULL; + } + +/* debugging.. + lwnotice("Bounds poly: %s", + lwgeom_to_ewkt(GEOS2LWGEOM(g1_bounds, 0), + PARSER_CHECK_NONE)); + lwnotice("Line: %s", + lwgeom_to_ewkt(GEOS2LWGEOM(g2, 0), + PARSER_CHECK_NONE)); + + lwnotice("Noded bounds: %s", + lwgeom_to_ewkt(GEOS2LWGEOM(vgeoms[0], 0), + PARSER_CHECK_NONE)); +*/ + + polygons = GEOSPolygonize(vgeoms, 1); + if ( NULL == polygons ) + { + GEOSGeom_destroy(g1); + GEOSGeom_destroy(g2); + GEOSGeom_destroy(g1_bounds); + GEOSGeom_destroy((GEOSGeometry*)vgeoms[0]); + lwerror("GEOSPolygonize: %s", lwgeom_geos_errmsg); + return NULL; + } + +#if PARANOIA_LEVEL > 0 + if ( GEOSGeometryTypeId(polygons) != COLLECTIONTYPE ) + { + GEOSGeom_destroy(g1); + GEOSGeom_destroy(g2); + GEOSGeom_destroy(g1_bounds); + GEOSGeom_destroy((GEOSGeometry*)vgeoms[0]); + GEOSGeom_destroy(polygons); + lwerror("Unexpected return from GEOSpolygonize"); + return 0; + } +#endif + + /* We should now have all polygons, just skip + * the ones which are in holes of the original + * geometries and return the rest in a collection + */ + n = GEOSGetNumGeometries(polygons); + out = lwcollection_construct(COLLECTIONTYPE, lwpoly_in->SRID, + NULL, 0, NULL); + /* Allocate space for all polys */ + out->geoms = lwalloc(sizeof(LWGEOM*)*n); + assert(0 == out->ngeoms); + for (i=0; igeoms[out->ngeoms++] = GEOS2LWGEOM(p, + TYPE_HASZ(lwpoly_in->type)); + } + + GEOSGeom_destroy(g1); + GEOSGeom_destroy(g2); + GEOSGeom_destroy(g1_bounds); + GEOSGeom_destroy((GEOSGeometry*)vgeoms[0]); + GEOSGeom_destroy(polygons); + + return (LWGEOM*)out; +} + +static LWGEOM* lwpoly_split(LWPOLY* lwpoly_in, LWGEOM* blade_in); static LWGEOM* -lwpoly_split(LWPOLY* lwgeom_in, LWGEOM* blade_in) +lwpoly_split(LWPOLY* lwpoly_in, LWGEOM* blade_in) { switch (TYPE_GETTYPE(blade_in->type)) { case LINETYPE: + return lwpoly_split_by_line(lwpoly_in, (LWLINE*)blade_in); default: lwerror("Splitting a Polygon by a %s is unsupported", lwgeom_typename(TYPE_GETTYPE(blade_in->type))); diff --git a/regress/split.sql b/regress/split.sql index 3dc13c11d..8985885af 100644 --- a/regress/split.sql +++ b/regress/split.sql @@ -30,3 +30,15 @@ select '7', st_asewkt(ST_SplitGeometry('SRID=10;LINESTRING(0 0, 10 0, 10 10, 0 1 select '8.1', st_asewkt(ST_SplitGeometry('SRID=10;LINESTRING(0 0, 10 0)', 'SRID=10;LINESTRING(5 0, 20 0)')); -- Split line by contained line (2) select '8.2', st_asewkt(ST_SplitGeometry('SRID=10;LINESTRING(0 0, 10 0)', 'SRID=10;LINESTRING(5 0, 8 0)')); + +-- Split exterior-only polygon by crossing line +select '20', st_asewkt(ST_SplitGeometry('SRID=12;POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))', 'SRID=12;LINESTRING(5 -5, 5 15)')); + +-- Split single-hole polygon by line crossing both exterior and hole +select '21', st_asewkt(ST_SplitGeometry('SRID=12;POLYGON((0 0, 10 0, 10 10, 0 10, 0 0),(2 2, 8 2, 8 8, 2 8, 2 2))', 'SRID=12;LINESTRING(5 -5, 5 15)')); + +-- Split single-hole polygon by line crossing only exterior +select '22', st_asewkt(ST_SplitGeometry('SRID=12;POLYGON((0 0, 10 0, 10 10, 0 10, 0 0),(5 2, 8 2, 8 8, 5 8, 5 2))', 'SRID=12;LINESTRING(2 -5, 2 15)')); + +-- Split double-hole polygon by line crossing exterior and both holes +select '23', st_asewkt(ST_SplitGeometry('SRID=12;POLYGON((0 0, 10 0, 10 10, 0 10, 0 0),(2 2, 8 2, 8 4, 2 4, 2 2),(2 6,8 6,8 8,2 8,2 6))', 'SRID=12;LINESTRING(5 -5, 5 15)')); diff --git a/regress/split_expected b/regress/split_expected index ea59b4196..c1d8fefad 100644 --- a/regress/split_expected +++ b/regress/split_expected @@ -10,3 +10,7 @@ ERROR: Operation on mixed SRID geometries 7|SRID=10;GEOMETRYCOLLECTION(LINESTRING(0 0,5 0),LINESTRING(5 0,10 0,10 10,5 10),LINESTRING(5 10,0 10,0 20,5 20),LINESTRING(5 20,10 20)) ERROR: Splitter line has linear intersection with input ERROR: Splitter line has linear intersection with input +20|SRID=12;GEOMETRYCOLLECTION(POLYGON((5 0,0 0,0 10,5 10,5 0)),POLYGON((5 10,10 10,10 0,5 0,5 10))) +21|SRID=12;GEOMETRYCOLLECTION(POLYGON((5 0,0 0,0 10,5 10,5 8,2 8,2 2,5 2,5 0)),POLYGON((5 10,10 10,10 0,5 0,5 2,8 2,8 8,5 8,5 10))) +22|SRID=12;GEOMETRYCOLLECTION(POLYGON((2 0,0 0,0 10,2 10,2 0)),POLYGON((2 10,10 10,10 0,2 0,2 10),(5 2,8 2,8 8,5 8,5 2))) +23|SRID=12;GEOMETRYCOLLECTION(POLYGON((5 0,0 0,0 10,5 10,5 8,2 8,2 6,5 6,5 4,2 4,2 2,5 2,5 0)),POLYGON((5 10,10 10,10 0,5 0,5 2,8 2,8 4,5 4,5 6,8 6,8 8,5 8,5 10))) -- 2.40.0