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

diff --git a/liblwgeom/lwgeom_geos.c b/liblwgeom/lwgeom_geos.c
index 4eebe51caa1600a5cd65a5f035ba0d0864290d69..491984efa1fca60b1f32ed92d1f34c18ed502c85 100644 (file)
--- 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; i<lwpoly->nrings; ++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; i<ngeoms; ++i)
                {
-                       GEOSGeometry* g = LWGEOM2GEOS(lwc->geoms[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 31f1eb990b1f189cce384f00e7ca6695608c90f9..a255d2b58dc2ae608237029666a803f0372ad198 100644 (file)
--- 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 71d340a3490d7426937e61863e1b93e94dd690db..3e0b5af878d146498e95829df28d67d5958d9cf7 100644 (file)
--- 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 3084c92d72415d401a079420b0a7a158d95fc54e..ac559fdeafc14f0698aba9bc1fcab0fa8ea68ac8 100644 (file)
--- 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 33e5890ac6201efc61d5e25da565dc5394384b4e..ee0828e8e19096974e8ac720a8c48e7d684d1e31 100644 (file)
--- 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 8c44e7a06fc787df9dc97397aafffb4abb5f0e15..1f863de34291fdea0b92c976473af5e6559c24e2 100644 (file)
--- 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 54b5dd27345790e741c950f73e738512fd9ff6b6..299e0b6eb5a6f2aac52c1bd1f38c20543685ca68 100644 (file)
--- 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))