Implement split-poly-by-line

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

diff --git a/postgis/lwgeom_geos_split.c b/postgis/lwgeom_geos_split.c
index 39698e82e1c923bb7429e85e46184005fae4ba27..f88986f9ce7640c3ab3809b8fc8b1d9d4e3e88e5 100644 (file)
--- 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; i<n; ++i)
+       {
+               GEOSGeometry* pos; /* point on surface */
+               const GEOSGeometry* p = GEOSGetGeometryN(polygons, i); 
+               int contains;
+
+               pos = GEOSPointOnSurface(p);
+               if ( ! pos ) {
+                       GEOSGeom_destroy(g1);
+                       GEOSGeom_destroy(g2);
+                       GEOSGeom_destroy(g1_bounds);
+                       GEOSGeom_destroy((GEOSGeometry*)vgeoms[0]);
+                       GEOSGeom_destroy(polygons);
+                       lwerror("GEOSPointOnSurface: %s", lwgeom_geos_errmsg);
+                       return NULL;
+               }
+
+               contains = GEOSContains(g1, pos);
+               if ( 2 == contains )
+               {
+                       GEOSGeom_destroy(g1);
+                       GEOSGeom_destroy(g2);
+                       GEOSGeom_destroy(g1_bounds);
+                       GEOSGeom_destroy((GEOSGeometry*)vgeoms[0]);
+                       GEOSGeom_destroy(polygons);
+                       GEOSGeom_destroy(pos);
+                       lwerror("GEOSContains: %s", lwgeom_geos_errmsg);
+                       return NULL;
+               }
+
+               GEOSGeom_destroy(pos);
+
+               if ( 0 == contains ) {
+                       /* Original geometry doesn't contain
+                        * a point in this ring, must be an hole
+                        */
+                       continue;
+               }
+
+               out->geoms[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 3dc13c11d3f5b4e6f149575b0564776233a7951a..8985885afd160eed07bc74c0b027de537951afe5 100644 (file)
--- 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 ea59b4196648cfa8b1b4f46f0bf30fb921bd296d..c1d8fefadc25c405b455e4c8bd4b543a42e63115 100644 (file)
--- 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)))