Do not drop polygon boundaries collapsed to points

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

diff --git a/postgis/lwgeom_geos_clean.c b/postgis/lwgeom_geos_clean.c
index 3953f932bed4bcad4ca5247ea1dc7f93a7d13ca6..f7796fbb1be9efef11b53ad59e322210ac42c5a2 100644 (file)
--- a/postgis/lwgeom_geos_clean.c
+++ b/postgis/lwgeom_geos_clean.c
@@ -391,8 +391,9 @@ LWGEOM_GEOS_makeValidPolygon(const GEOSGeometry* gin)
 {
        GEOSGeom gout;
        GEOSGeom geos_bound;
-       GEOSGeom geos_cut_edges, geos_area;
-       GEOSGeometry *vgeoms[2]; /* One for area, one for cut-edges */
+       GEOSGeom geos_cut_edges, geos_area, collapse_points;
+       GEOSGeometry *vgeoms[3]; /* One for area, one for cut-edges */
+       unsigned int nvgeoms=0;
 
        assert (GEOSGeomTypeId(gin) == GEOS_POLYGON ||
                GEOSGeomTypeId(gin) == GEOS_MULTIPOLYGON);
@@ -410,14 +411,72 @@ LWGEOM_GEOS_makeValidPolygon(const GEOSGeometry* gin)
 
        /* Use noded boundaries as initial "cut" edges */
 
+
        geos_cut_edges = LWGEOM_GEOS_nodeLines(geos_bound);
-       GEOSGeom_destroy(geos_bound); 
        if ( NULL == geos_cut_edges )
        {
+               GEOSGeom_destroy(geos_bound); 
                lwnotice("LWGEOM_GEOS_nodeLines(): %s", lwgeom_geos_errmsg);
                return NULL;
        }
 
+       /* NOTE: the noding process may drop lines collapsing to points.
+        *       We want to retrive any of those */
+       {
+               GEOSGeometry* pi;
+               GEOSGeometry* po;
+               
+               pi = GEOSGeom_extractUniquePoints(geos_bound);
+               if ( NULL == pi ) {
+                       GEOSGeom_destroy(geos_bound); 
+                       lwnotice("GEOSGeom_extractUniquePoints(): %s",
+                               lwgeom_geos_errmsg);
+                       return NULL;
+               }
+
+       POSTGIS_DEBUGF(3,
+                      "Boundaries input points %s",
+                      lwgeom_to_ewkt(GEOS2LWGEOM(pi, 0),
+                                     PARSER_CHECK_NONE));
+
+               po = GEOSGeom_extractUniquePoints(geos_cut_edges);
+               if ( NULL == po ) {
+                       GEOSGeom_destroy(geos_bound); 
+                       GEOSGeom_destroy(pi); 
+                       lwnotice("GEOSGeom_extractUniquePoints(): %s",
+                               lwgeom_geos_errmsg);
+                       return NULL;
+               }
+
+       POSTGIS_DEBUGF(3,
+                      "Boundaries output points %s",
+                      lwgeom_to_ewkt(GEOS2LWGEOM(po, 0),
+                                     PARSER_CHECK_NONE));
+
+               collapse_points = GEOSDifference(pi, po);
+               if ( NULL == collapse_points ) {
+                       GEOSGeom_destroy(geos_bound); 
+                       GEOSGeom_destroy(pi); 
+                       GEOSGeom_destroy(po); 
+                       lwnotice("GEOSDifference(): %s", lwgeom_geos_errmsg);
+                       return NULL;
+               }
+
+       POSTGIS_DEBUGF(3,
+                      "Collapse points: %s",
+                      lwgeom_to_ewkt(GEOS2LWGEOM(collapse_points, 0),
+                                     PARSER_CHECK_NONE));
+
+               GEOSGeom_destroy(pi);
+               GEOSGeom_destroy(po);
+       }
+       GEOSGeom_destroy(geos_bound); 
+
+       POSTGIS_DEBUGF(3,
+                      "Noded Boundaries: %s",
+                      lwgeom_to_ewkt(GEOS2LWGEOM(geos_cut_edges, 0),
+                                     PARSER_CHECK_NONE));
+
        /* And use an empty geometry as initial "area" */
        geos_area = GEOSGeom_createEmptyPolygon();
        if ( ! geos_area )
@@ -524,34 +583,45 @@ LWGEOM_GEOS_makeValidPolygon(const GEOSGeometry* gin)
                geos_cut_edges = new_cut_edges;
        }
 
-       if ( GEOSisEmpty(geos_area) )
+       if ( ! GEOSisEmpty(geos_area) ) {
+               vgeoms[nvgeoms++] = geos_area;
+       } else {
+               GEOSGeom_destroy(geos_area);
+       }
+
+       if ( ! GEOSisEmpty(geos_cut_edges) ) {
+               vgeoms[nvgeoms++] = geos_cut_edges;
+       } else {
+               GEOSGeom_destroy(geos_cut_edges);
+       }
+
+       if ( ! GEOSisEmpty(collapse_points) ) {
+               vgeoms[nvgeoms++] = collapse_points;
+       } else {
+               GEOSGeom_destroy(collapse_points);
+       }
+
+       if ( 1 == nvgeoms )
        {
                /* Return cut edges */
-               GEOSGeom_destroy(geos_area);
-               return geos_cut_edges;
+               gout = vgeoms[0];
        }
-       else if ( ! GEOSisEmpty(geos_cut_edges) )
+       else 
        {
                /* Collect areas and lines (if any line) */
-               vgeoms[0] = geos_area;
-               vgeoms[1] = geos_cut_edges;
-               gout = GEOSGeom_createCollection(GEOS_GEOMETRYCOLLECTION, vgeoms, 2);
+               gout = GEOSGeom_createCollection(GEOS_GEOMETRYCOLLECTION, vgeoms, nvgeoms);
                if ( ! gout )   /* an exception again */
                {
                        /* cleanup and throw */
                        lwnotice("GEOSGeom_createCollection() threw an error: %s",
                                 lwgeom_geos_errmsg);
+                       /* TODO: cleanup! */
                        return NULL;
                }
-
-               return gout;
-       }
-       else
-       {
-               GEOSGeom_destroy(geos_cut_edges);
-               return geos_area;
        }
 
+       return gout;
+
 }
 
 static GEOSGeometry*

diff --git a/regress/clean.sql b/regress/clean.sql
index 2f1905ec8d389ee3d5cce39ff9bb1a536fb25909..218e023756061e4cfd4a5a3fac9bacc86fac1fcd 100644 (file)
--- a/regress/clean.sql
+++ b/regress/clean.sql
@@ -28,6 +28,7 @@ RT    17.1    010300000003000000050000000000000000004E4000000000000014C000000000000054
 PG     1       0103000000010000000100000000000000000000000000000000000000      POINT(0 0)
 PG     2       LINESTRING(0 0, 0 0)    POINT(0 0)
 PG     3       MULTILINESTRING((0 0, 10 0),(20 20, 20 20))     GEOMETRYCOLLECTION(LINESTRING(0 0, 10 0),POINT(20 20))
+PG     4       MULTIPOLYGON(((5 3, 7 4, 9 5, 11 6, 13 7, 5 3)),((14 14, 14 14, 14 14, 14 14))) GEOMETRYCOLLECTION(MULTILINESTRING((5 3,7 4),(7 4,9 5),(9 5,11 6),(11 6,13 7)),POINT(14 14))
 \.
 
 -- PG.1 : polygon with single ring with single point in it

diff --git a/regress/clean_expected b/regress/clean_expected
index 3cfbb9f935ef85199377f685b15e62022c0e7e6b..5f6663b4998c689e67d36fabd4d363cc62e1fb40 100644 (file)
--- a/regress/clean_expected
+++ b/regress/clean_expected
@@ -24,6 +24,7 @@ RT|17.1|t|t|f
 PG|1|t|t|f
 PG|2|t|t|f
 PG|3|t|t|f
+PG|4|t|t|f
 SRID1|3
 SRID2|3
 SRID3|3