From ba2c3d1b60d69a7f74bbb57296e31347dc29f89a Mon Sep 17 00:00:00 2001
From: Sandro Santilli <strk@keybit.net>
Date: Mon, 3 May 2010 11:31:18 +0000
Subject: [PATCH] Do not drop polygon boundaries collapsed to points

git-svn-id: http://svn.osgeo.org/postgis/trunk@5599 b70326c6-7e19-0410-871a-916f4a2858ee
---
 postgis/lwgeom_geos_clean.c | 104 ++++++++++++++++++++++++++++++------
 regress/clean.sql           |   1 +
 regress/clean_expected      |   1 +
 3 files changed, 89 insertions(+), 17 deletions(-)

diff --git a/postgis/lwgeom_geos_clean.c b/postgis/lwgeom_geos_clean.c
index 3953f932b..f7796fbb1 100644
--- 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 2f1905ec8..218e02375 100644
--- 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 3cfbb9f93..5f6663b49 100644
--- 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
-- 
2.40.0