From: Sandro Santilli <strk@keybit.net>
Date: Thu, 25 Feb 2010 13:15:55 +0000 (+0000)
Subject: Recursively try to build area with portions of the original boundary not on the bound... 
X-Git-Tag: 2.0.0alpha1~3176
X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=500202378487f3d6ba6c2b1bb80decef9fe43fb5;p=postgis

Recursively try to build area with portions of the original boundary not on the boundary of newly constructed area and symdifference the new area (if any) with the final polygon being built. This behaviour gives more chance to get a pure areal (not mixed) output still not missing vertices.


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

diff --git a/postgis/lwgeom_geos_clean.c b/postgis/lwgeom_geos_clean.c
index c98b64f54..b23617b02 100644
--- a/postgis/lwgeom_geos_clean.c
+++ b/postgis/lwgeom_geos_clean.c
@@ -478,6 +478,7 @@ LWGEOM_GEOS_makeValidPolygon(const GEOSGeometry* gin)
 		         lwgeom_to_ewkt(GEOS2LWGEOM(geos_area, 0),
 		                        PARSER_CHECK_NONE),
 		         loggederror);
+		GEOSGeom_destroy(geos_area);
 		return geos_bound_noded;
 	}
 
@@ -524,12 +525,108 @@ LWGEOM_GEOS_makeValidPolygon(const GEOSGeometry* gin)
 		return geos_area;
 	}
 
+	/*
+	 * See if an area can be build with the remaining edges
+	 * and if it can, symdifference with the original area.
+	 * Iterate this until no more polygons can be created
+	 * with left-over edges.
+	 */
+	while (GEOSGetNumGeometries(geos_cut_edges))
+	{
+		GEOSGeometry* new_area=0;
+		GEOSGeometry* new_area_bound=0;
+		GEOSGeometry* symdif=0;
+		GEOSGeometry* new_cut_edges=0;
+
+		/*
+		 * ASSUMPTION: cut_edges should already be fully noded
+		 */
+
+		new_area = LWGEOM_GEOS_buildArea(geos_cut_edges);
+		if ( ! new_area ) { /* must be an exception */
+			GEOSGeom_destroy(geos_cut_edges);
+			GEOSGeom_destroy(geos_area);
+			lwnotice("LWGEOM_GEOS_buildArea() threw an error: %s",
+				 loggederror);
+			return NULL;
+		}
+
+		if ( GEOSisEmpty(new_area) ) {
+			/* no more rings can be build with thes edges */
+			GEOSGeom_destroy(new_area);
+			break;
+		}
+
+		/*
+		 * We succeeded in building a ring !
+		 */
+
+
+		/*
+		 * Save the new ring boundaries first (to compute
+		 * further cut edges later)
+		 */
+		new_area_bound = GEOSBoundary(new_area);
+		if ( ! new_area_bound ) { 
+			/* We did check for empty area already so
+			 * this must be some other error */
+			lwnotice("GEOSBoundary('%s') threw an error: %s",
+				 lwgeom_to_ewkt(GEOS2LWGEOM(new_area, 0),
+						PARSER_CHECK_NONE),
+				 loggederror);
+			GEOSGeom_destroy(new_area);
+			GEOSGeom_destroy(geos_area);
+			return NULL;
+		}
+
+		/*
+		 * Now symdif new and old area 
+		 */
+		symdif = GEOSSymDifference(geos_area, new_area);
+		if ( ! symdif ) { /* must be an exception */
+			GEOSGeom_destroy(geos_cut_edges);
+			GEOSGeom_destroy(new_area);
+			GEOSGeom_destroy(new_area_bound);
+			GEOSGeom_destroy(geos_area);
+			lwnotice("GEOSSymDifference() threw an error: %s",
+				 loggederror);
+			return NULL;
+		}
+
+		GEOSGeom_destroy(geos_area);
+		GEOSGeom_destroy(new_area);
+		geos_area = symdif;
+		symdif = 0;
+
+		/*
+		 * Now let's re-set geos_cut_edges with what's left
+		 * from the original boundary.
+		 * ASSUMPTION: only the previous cut-edges can be
+		 *             left, so we don't need to reconsider
+		 *             the whole original boundaries
+		 */
+
+		new_cut_edges = GEOSDifference(geos_cut_edges, new_area_bound);
+		GEOSGeom_destroy(new_area_bound);
+		if ( ! new_cut_edges )   /* an exception ? */
+		{
+			/* cleanup and throw */
+			GEOSGeom_destroy(geos_cut_edges);
+			GEOSGeom_destroy(geos_area);
+			lwnotice("GEOSDifference() threw an error: %s",
+				 loggederror);
+			return NULL;
+		}
+		GEOSGeom_destroy(geos_cut_edges);
+		geos_cut_edges = new_cut_edges;
+	}
+
 	/*
 	 * Drop cut-edge segments for which all
 	 * vertices are already vertices of geos_bound_noded
 	 * which fall completely in the area
 	 */
-	{
+	if ( 0 ) {
 
 		GEOSGeometry** cutlines;
 		GEOSGeometry* geos_bound_noded_points;
@@ -795,6 +892,8 @@ lwgeom_make_valid(LWGEOM* lwgeom_in)
 		lwgeom_out = lwgeom_as_multi(lwgeom_out);
 	}
 
+	GEOSGeom_destroy(geosgeom);
+
 	return lwgeom_out;
 }