From b58ad4f2bb71fd02e7f8b92c80792cec58579d6c Mon Sep 17 00:00:00 2001 From: Sandro Santilli Date: Wed, 1 Oct 2014 12:35:58 +0000 Subject: [PATCH] Add support for auto-fix geom structure for GEOS conversions Fixes ST_ClipByBox2D run with non-closed polygon rings in input (#2945) git-svn-id: http://svn.osgeo.org/postgis/trunk@13018 b70326c6-7e19-0410-871a-916f4a2858ee --- liblwgeom/lwgeom_geos.c | 75 ++++++++++++++++++++++------------ liblwgeom/lwgeom_geos.h | 2 +- liblwgeom/lwgeom_geos_clean.c | 5 ++- postgis/lwgeom_geos.c | 8 ++-- postgis/lwgeom_geos_prepared.c | 2 +- regress/clipbybox2d.sql | 5 +++ regress/clipbybox2d_expected | 1 + 7 files changed, 65 insertions(+), 33 deletions(-) diff --git a/liblwgeom/lwgeom_geos.c b/liblwgeom/lwgeom_geos.c index 4eebe51ca..491984efa 100644 --- a/liblwgeom/lwgeom_geos.c +++ b/liblwgeom/lwgeom_geos.c @@ -234,11 +234,38 @@ ptarray_to_GEOSCoordSeq(const POINTARRAY *pa) return sq; } +static GEOSGeometry * +ptarray_to_GEOSLinearRing(const POINTARRAY *pa, int autofix) +{ + GEOSCoordSeq sq; + GEOSGeom g; + POINTARRAY *npa = 0; + if ( autofix ) + { + /* check ring for being closed and fix if not */ + if ( ! ptarray_is_closed_2d(pa) ) { + npa = ptarray_addPoint(pa, getPoint_internal(pa, 0), + FLAGS_NDIMS(pa->flags), pa->npoints); + pa = npa; + } + /* TODO: check ring for having at least 4 vertices */ +#if 0 + while ( pa->npoints < 4 ) { + npa = ptarray_addPoint(npa, getPoint_internal(pa, 0), + FLAGS_NDIMS(pa->flags), pa->npoints); + } +#endif + } + sq = ptarray_to_GEOSCoordSeq(pa); + if ( npa ) ptarray_free(npa); + g = GEOSGeom_createLinearRing(sq); + return g; +} GEOSGeometry * -LWGEOM2GEOS(const LWGEOM *lwgeom) +LWGEOM2GEOS(const LWGEOM *lwgeom, int autofix) { GEOSCoordSeq sq; GEOSGeom g, shell; @@ -299,6 +326,7 @@ LWGEOM2GEOS(const LWGEOM *lwgeom) break; case LINETYPE: lwl = (LWLINE *)lwgeom; + /* TODO: if (autofix) */ if ( lwl->points->npoints == 1 ) { /* Duplicate point, to make geos-friendly */ lwl->points = ptarray_addPoint(lwl->points, @@ -330,9 +358,7 @@ LWGEOM2GEOS(const LWGEOM *lwgeom) } else { - sq = ptarray_to_GEOSCoordSeq(lwpoly->rings[0]); - /* TODO: check ring for being closed and fix if not */ - shell = GEOSGeom_createLinearRing(sq); + shell = ptarray_to_GEOSLinearRing(lwpoly->rings[0], autofix); if ( ! shell ) return NULL; /*lwerror("LWGEOM2GEOS: exception during polygon shell conversion"); */ ngeoms = lwpoly->nrings-1; @@ -341,8 +367,7 @@ LWGEOM2GEOS(const LWGEOM *lwgeom) for (i=1; inrings; ++i) { - sq = ptarray_to_GEOSCoordSeq(lwpoly->rings[i]); - geoms[i-1] = GEOSGeom_createLinearRing(sq); + geoms[i-1] = ptarray_to_GEOSLinearRing(lwpoly->rings[i], autofix); if ( ! geoms[i-1] ) { --i; @@ -379,7 +404,7 @@ LWGEOM2GEOS(const LWGEOM *lwgeom) for (i=0; igeoms[i]); + GEOSGeometry* g = LWGEOM2GEOS(lwc->geoms[i], 0); if ( ! g ) { while (i) GEOSGeom_destroy(geoms[--i]); @@ -429,7 +454,7 @@ lwgeom_normalize(const LWGEOM *geom1) initGEOS(lwnotice, lwgeom_geos_error); - g1 = LWGEOM2GEOS(geom1); + g1 = LWGEOM2GEOS(geom1, 0); if ( 0 == g1 ) /* exception thrown at construction */ { lwerror("First argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg); @@ -482,14 +507,14 @@ lwgeom_intersection(const LWGEOM *geom1, const LWGEOM *geom2) LWDEBUG(3, "intersection() START"); - g1 = LWGEOM2GEOS(geom1); + g1 = LWGEOM2GEOS(geom1, 0); if ( 0 == g1 ) /* exception thrown at construction */ { lwerror("First argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg); return NULL ; } - g2 = LWGEOM2GEOS(geom2); + g2 = LWGEOM2GEOS(geom2, 0); if ( 0 == g2 ) /* exception thrown at construction */ { lwerror("Second argument geometry could not be converted to GEOS."); @@ -563,14 +588,14 @@ lwgeom_difference(const LWGEOM *geom1, const LWGEOM *geom2) initGEOS(lwnotice, lwgeom_geos_error); - g1 = LWGEOM2GEOS(geom1); + g1 = LWGEOM2GEOS(geom1, 0); if ( 0 == g1 ) /* exception thrown at construction */ { lwerror("First argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg); return NULL; } - g2 = LWGEOM2GEOS(geom2); + g2 = LWGEOM2GEOS(geom2, 0); if ( 0 == g2 ) /* exception thrown at construction */ { GEOSGeom_destroy(g1); @@ -637,7 +662,7 @@ lwgeom_symdifference(const LWGEOM* geom1, const LWGEOM* geom2) initGEOS(lwnotice, lwgeom_geos_error); - g1 = LWGEOM2GEOS(geom1); + g1 = LWGEOM2GEOS(geom1, 0); if ( 0 == g1 ) /* exception thrown at construction */ { @@ -645,7 +670,7 @@ lwgeom_symdifference(const LWGEOM* geom1, const LWGEOM* geom2) return NULL; } - g2 = LWGEOM2GEOS(geom2); + g2 = LWGEOM2GEOS(geom2, 0); if ( 0 == g2 ) /* exception thrown at construction */ { @@ -713,7 +738,7 @@ lwgeom_union(const LWGEOM *geom1, const LWGEOM *geom2) initGEOS(lwnotice, lwgeom_geos_error); - g1 = LWGEOM2GEOS(geom1); + g1 = LWGEOM2GEOS(geom1, 0); if ( 0 == g1 ) /* exception thrown at construction */ { @@ -721,7 +746,7 @@ lwgeom_union(const LWGEOM *geom1, const LWGEOM *geom2) return NULL; } - g2 = LWGEOM2GEOS(geom2); + g2 = LWGEOM2GEOS(geom2, 0); if ( 0 == g2 ) /* exception thrown at construction */ { @@ -787,7 +812,7 @@ lwgeom_clip_by_rect(const LWGEOM *geom1, double x0, double y0, double x1, double LWDEBUG(3, "clip_by_rect() START"); - g1 = LWGEOM2GEOS(geom1); + g1 = LWGEOM2GEOS(geom1, 1); /* auto-fix structure */ if ( 0 == g1 ) /* exception thrown at construction */ { lwerror("First argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg); @@ -1123,7 +1148,7 @@ lwgeom_buildarea(const LWGEOM *geom) initGEOS(lwnotice, lwgeom_geos_error); - geos_in = LWGEOM2GEOS(geom); + geos_in = LWGEOM2GEOS(geom, 0); if ( 0 == geos_in ) /* exception thrown at construction */ { @@ -1172,7 +1197,7 @@ lwgeom_geos_noop(const LWGEOM* geom_in) int is3d = FLAGS_GET_Z(geom_in->flags); initGEOS(lwnotice, lwgeom_geos_error); - geosgeom = LWGEOM2GEOS(geom_in); + geosgeom = LWGEOM2GEOS(geom_in, 0); if ( ! geosgeom ) { lwerror("Geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg); @@ -1210,14 +1235,14 @@ lwgeom_snap(const LWGEOM* geom1, const LWGEOM* geom2, double tolerance) initGEOS(lwnotice, lwgeom_geos_error); - g1 = (GEOSGeometry *)LWGEOM2GEOS(geom1); + g1 = (GEOSGeometry *)LWGEOM2GEOS(geom1, 0); if ( 0 == g1 ) /* exception thrown at construction */ { lwerror("First argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg); return NULL; } - g2 = (GEOSGeometry *)LWGEOM2GEOS(geom2); + g2 = (GEOSGeometry *)LWGEOM2GEOS(geom2, 0); if ( 0 == g2 ) /* exception thrown at construction */ { lwerror("Second argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg); @@ -1273,14 +1298,14 @@ lwgeom_sharedpaths(const LWGEOM* geom1, const LWGEOM* geom2) initGEOS(lwnotice, lwgeom_geos_error); - g1 = (GEOSGeometry *)LWGEOM2GEOS(geom1); + g1 = (GEOSGeometry *)LWGEOM2GEOS(geom1, 0); if ( 0 == g1 ) /* exception thrown at construction */ { lwerror("First argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg); return NULL; } - g2 = (GEOSGeometry *)LWGEOM2GEOS(geom2); + g2 = (GEOSGeometry *)LWGEOM2GEOS(geom2, 0); if ( 0 == g2 ) /* exception thrown at construction */ { lwerror("Second argument geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg); @@ -1325,7 +1350,7 @@ lwgeom_offsetcurve(const LWLINE *lwline, double size, int quadsegs, int joinStyl initGEOS(lwnotice, lwgeom_geos_error); - g1 = (GEOSGeometry *)LWGEOM2GEOS(lwgeom_in); + g1 = (GEOSGeometry *)LWGEOM2GEOS(lwgeom_in, 0); if ( ! g1 ) { lwerror("lwgeom_offsetcurve: Geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg); @@ -1446,7 +1471,7 @@ LWGEOM* lwgeom_delaunay_triangulation(const LWGEOM *lwgeom_in, double tolerance, initGEOS(lwnotice, lwgeom_geos_error); - g1 = (GEOSGeometry *)LWGEOM2GEOS(lwgeom_in); + g1 = (GEOSGeometry *)LWGEOM2GEOS(lwgeom_in, 0); if ( ! g1 ) { lwerror("lwgeom_delaunay_triangulation: Geometry could not be converted to GEOS: %s", lwgeom_geos_errmsg); diff --git a/liblwgeom/lwgeom_geos.h b/liblwgeom/lwgeom_geos.h index 31f1eb990..a255d2b58 100644 --- a/liblwgeom/lwgeom_geos.h +++ b/liblwgeom/lwgeom_geos.h @@ -24,7 +24,7 @@ ** Public prototypes for GEOS utility functions. */ LWGEOM *GEOS2LWGEOM(const GEOSGeometry *geom, char want3d); -GEOSGeometry * LWGEOM2GEOS(const LWGEOM *g); +GEOSGeometry * LWGEOM2GEOS(const LWGEOM *g, int autofix); GEOSGeometry * LWGEOM_GEOS_buildArea(const GEOSGeometry* geom_in); diff --git a/liblwgeom/lwgeom_geos_clean.c b/liblwgeom/lwgeom_geos_clean.c index 71d340a34..3e0b5af87 100644 --- a/liblwgeom/lwgeom_geos_clean.c +++ b/liblwgeom/lwgeom_geos_clean.c @@ -1009,7 +1009,7 @@ lwgeom_make_valid(LWGEOM* lwgeom_in) initGEOS(lwgeom_geos_error, lwgeom_geos_error); lwgeom_out = lwgeom_in; - geosgeom = LWGEOM2GEOS(lwgeom_out); + geosgeom = LWGEOM2GEOS(lwgeom_out, 0); if ( ! geosgeom ) { LWDEBUGF(4, @@ -1025,7 +1025,8 @@ lwgeom_make_valid(LWGEOM* lwgeom_in) } /* try again as we did cleanup now */ - geosgeom = LWGEOM2GEOS(lwgeom_out); + /* TODO: invoke LWGEOM2GEOS directly with autoclean ? */ + geosgeom = LWGEOM2GEOS(lwgeom_out, 0); if ( ! geosgeom ) { lwerror("Couldn't convert POSTGIS geom to GEOS: %s", diff --git a/postgis/lwgeom_geos.c b/postgis/lwgeom_geos.c index 3084c92d7..ac559fdea 100644 --- a/postgis/lwgeom_geos.c +++ b/postgis/lwgeom_geos.c @@ -586,7 +586,7 @@ Datum pgis_union_geometry_array(PG_FUNCTION_ARGS) geoms = repalloc( geoms, sizeof(GEOSGeom) * geoms_size ); } /* This builds a LWPOLY on top of the serialized form */ - g = LWGEOM2GEOS(lwpoly_as_lwgeom(lwmpoly->geoms[j])); + g = LWGEOM2GEOS(lwpoly_as_lwgeom(lwmpoly->geoms[j], 0)); if ( 0 == g ) /* exception thrown at construction */ { /* TODO: cleanup all GEOS memory */ @@ -928,7 +928,7 @@ Datum boundary(PG_FUNCTION_ARGS) initGEOS(lwnotice, lwgeom_geos_error); - g1 = LWGEOM2GEOS(lwgeom); + g1 = LWGEOM2GEOS(lwgeom, 0); lwgeom_free(lwgeom); if ( 0 == g1 ) /* exception thrown at construction */ @@ -1785,7 +1785,7 @@ Datum isvalid(PG_FUNCTION_ARGS) { lwerror("unable to deserialize input"); } - g1 = LWGEOM2GEOS(lwgeom); + g1 = LWGEOM2GEOS(lwgeom, 0); lwgeom_free(lwgeom); if ( ! g1 ) @@ -3264,7 +3264,7 @@ POSTGIS2GEOS(GSERIALIZED *pglwgeom) lwerror("POSTGIS2GEOS: unable to deserialize input"); return NULL; } - ret = LWGEOM2GEOS(lwgeom); + ret = LWGEOM2GEOS(lwgeom, 0); lwgeom_free(lwgeom); if ( ! ret ) { diff --git a/postgis/lwgeom_geos_prepared.c b/postgis/lwgeom_geos_prepared.c index 33e5890ac..ee0828e8e 100644 --- a/postgis/lwgeom_geos_prepared.c +++ b/postgis/lwgeom_geos_prepared.c @@ -318,7 +318,7 @@ PrepGeomCacheBuilder(const LWGEOM *lwgeom, GeomCache *cache) return LW_FAILURE; } - prepcache->geom = LWGEOM2GEOS( lwgeom ); + prepcache->geom = LWGEOM2GEOS( lwgeom , 0); if ( ! prepcache->geom ) return LW_FAILURE; prepcache->prepared_geom = GEOSPrepare( prepcache->geom ); if ( ! prepcache->prepared_geom ) return LW_FAILURE; diff --git a/regress/clipbybox2d.sql b/regress/clipbybox2d.sql index 8c44e7a06..1f863de34 100644 --- a/regress/clipbybox2d.sql +++ b/regress/clipbybox2d.sql @@ -10,4 +10,9 @@ SELECT ST_AsText(ST_ClipByBox2d('MULTIPOINT(-1 -1, 0 0, 2 2)'::geometry,ST_MakeE SELECT ST_AsText(ST_ClipByBox2d('POLYGON((0 0, 10 10, 0 10, 10 0, 0 0))', ST_MakeEnvelope(2,2,8,8))); -- Invalid polygon (lineal self-intersection) -- ST_Intersection returns a collection SELECT ST_AsText(ST_ClipByBox2d('POLYGON((0 0,5 4,5 6,0 10,10 10,5 6,5 4,10 0,0 0))', ST_MakeEnvelope(2,2,10,5))); +-- Invalid polygon (non-closed ring) -- ST_Intersection raises an exception +-- The HEXWKB represents 'POLYGON((0 0, 10 0, 10 10, 0 10))' +SELECT ST_AsText(ST_ClipByBox2d( + '0103000000010000000400000000000000000000000000000000000000000000000000244000000000000000000000000000002440000000000000244000000000000000000000000000002440' + , ST_MakeEnvelope(2,2,5,5))); diff --git a/regress/clipbybox2d_expected b/regress/clipbybox2d_expected index 54b5dd273..299e0b6eb 100644 --- a/regress/clipbybox2d_expected +++ b/regress/clipbybox2d_expected @@ -4,3 +4,4 @@ BOX(2 2,8 8) POINT(2 2) POLYGON((2 2,8 2,2 8,8 8,2 2)) POLYGON((2.5 2,5 4,5 5,5 4,7.5 2,2.5 2)) +POLYGON((2 2,2 5,5 5,5 2,2 2)) -- 2.50.1