From: Sandro Santilli <strk@keybit.net>
Date: Wed, 9 May 2012 12:08:16 +0000 (+0000)
Subject: Significatively speedup BuildArea with complex input (#1806)
X-Git-Tag: 2.0.1~85
X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=f81decb8a14976de732ea833c2c26937202b1843;p=postgis

Significatively speedup BuildArea with complex input (#1806)

Affects ST_BuildArea, ST_MakeValid and ST_GetFaceGeometry.

Replaces the iterated SymDifference used since 2005 with a more
scalable algorithm. The new algorithm removes from the polygonized
result all polygons whose rings are known to be already represented
by holes or shells of other polygons and finally runs a single
overlay operation (unary union).

With the case attached to ticket #1806, ST_BuildArea completes
within 12 seconds using the new code while it takes 27 _minutes_
with the old. Both versions return the same result (according to
ST_Equals).

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

diff --git a/NEWS b/NEWS
index 109e3d676..22844727d 100644
--- a/NEWS
+++ b/NEWS
@@ -17,12 +17,13 @@ PostGIS 2.0.1
   - #1789, fix false edge-node crossing report in ValidateTopology.
   - #1790, fix toTopoGeom handling of duplicated primitives.
   - #1791, fix ST_Azimuth with very close but distinct points.
-  - #1802, improved function interruptibility.
 
 * Enhancements *
 
   - More detailed exception message from topology editing functions.
   - #1786, improved build dependencies 
+  - #1802, improved function interruptibility.
+  - #1806, speedup of ST_BuildArea, ST_MakeValid and ST_GetFaceGeometry.
 
 PostGIS 2.0.0
 2012/04/03
diff --git a/liblwgeom/cunit/cu_buildarea.c b/liblwgeom/cunit/cu_buildarea.c
index eefbf98fe..a8a090b2d 100644
--- a/liblwgeom/cunit/cu_buildarea.c
+++ b/liblwgeom/cunit/cu_buildarea.c
@@ -300,7 +300,7 @@ static void buildarea7(void)
 	CU_ASSERT(gin != NULL);
 
 	gexp = lwgeom_from_wkt(
-"MULTIPOLYGON(((80 0,80 70,110 70,110 0,80 0),(90 60,90 50,100 50,100 60,90 60)),((20 20,20 30,20 50,30 50,30 30,30 20,20 20)),((0 0,0 70,70 70,70 0,0 0),(50 20,60 20,60 40,60 60,50 60,50 40,50 20),(10 10,40 10,40 60,10 60,10 10)))",
+"MULTIPOLYGON(((80 0,80 70,110 70,110 0,80 0),(90 60,90 50,100 50,100 60,90 60)),((20 20,20 30,20 50,30 50,30 30,30 20,20 20)),((0 0,0 70,70 70,70 0,0 0),(10 10,40 10,40 60,10 60,10 10),(50 20,60 20,60 40,60 60,50 60,50 40,50 20)))",
 		LW_PARSER_CHECK_NONE);
 	CU_ASSERT(gexp != NULL);
 
diff --git a/liblwgeom/lwgeom_geos.c b/liblwgeom/lwgeom_geos.c
index 2ea4efb33..d968f0d4d 100644
--- a/liblwgeom/lwgeom_geos.c
+++ b/liblwgeom/lwgeom_geos.c
@@ -721,140 +721,209 @@ lwgeom_union(const LWGEOM *geom1, const LWGEOM *geom2)
 	return result;
 }
 
-static int
-compare_by_numpoints(const void* g1, const void* g2)
+/* ------------ BuildArea stuff ---------------------------------------------------------------------{ */
+
+typedef struct Face_t {
+  const GEOSGeometry* geom;
+  GEOSGeometry* env;
+  double envarea;
+  struct Face_t* parent; /* if this face is an hole of another one, or NULL */
+} Face;
+
+static Face* newFace(const GEOSGeometry* g);
+static void delFace(Face* f);
+static unsigned int countParens(const Face* f);
+static void findFaceHoles(Face** faces, int nfaces);
+
+static Face*
+newFace(const GEOSGeometry* g)
 {
-  int n1 = GEOSGetNumCoordinates(*(const GEOSGeometry**)g1);
-  int n2 = GEOSGetNumCoordinates(*(const GEOSGeometry**)g2);
-  if ( n1 < n2 ) return -1;
-  if ( n1 > n2 ) return 1;
-  return 0;
+  Face* f = lwalloc(sizeof(Face));
+  f->geom = g;
+  f->env = GEOSEnvelope(f->geom);
+  GEOSArea(f->env, &f->envarea);
+  f->parent = NULL;
+  /* lwnotice("Built Face with area %g and %d holes", f->envarea, GEOSGetNumInteriorRings(f->geom)); */
+  return f;
 }
 
-GEOSGeometry*
-LWGEOM_GEOS_buildArea(const GEOSGeometry* geom_in)
+static unsigned int
+countParens(const Face* f)
 {
-	GEOSGeometry *tmp;
-	GEOSGeometry *geos_result, *shp;
-	GEOSGeometry const *vgeoms[1];
-	uint32_t i, ngeoms;
-	int srid = GEOSGetSRID(geom_in);
-	const GEOSGeometry ** geoms;
-
-	vgeoms[0] = geom_in;
-	geos_result = GEOSPolygonize(vgeoms, 1);
+  unsigned int pcount = 0;
+  while ( f->parent ) {
+    ++pcount;
+    f = f->parent;
+  }
+  return pcount;
+}
 
-	LWDEBUGF(3, "GEOSpolygonize returned @ %p", geos_result);
+/* Destroy the face and release memory associated with it */
+static void
+delFace(Face* f)
+{
+  GEOSGeom_destroy(f->env);
+  lwfree(f);
+}
 
-	/* Null return from GEOSpolygonize (an exception) */
-	if ( ! geos_result ) return 0;
 
-	/*
-	 * We should now have a collection
-	 */
-#if PARANOIA_LEVEL > 0
-	if ( GEOSGeometryTypeId(geos_result) != COLLECTIONTYPE )
-	{
-		GEOSGeom_destroy(geos_result);
-		lwerror("Unexpected return from GEOSpolygonize");
-		return 0;
-	}
-#endif
-
-	ngeoms = GEOSGetNumGeometries(geos_result);
+static int
+compare_by_envarea(const void* g1, const void* g2)
+{
+  Face* f1 = *(Face**)g1;
+  Face* f2 = *(Face**)g2;
+  double n1 = f1->envarea;
+  double n2 = f2->envarea;
 
-	LWDEBUGF(3, "GEOSpolygonize: ngeoms in polygonize output: %d", ngeoms);
-	LWDEBUGF(3, "GEOSpolygonize: polygonized:%s",
-	               lwgeom_to_ewkt(GEOS2LWGEOM(geos_result, 0)));
+  if ( n1 < n2 ) return 1;
+  if ( n1 > n2 ) return -1;
+  return 0;
+}
 
-	/*
-	 * No geometries in collection, early out
-	 */
-	if ( ngeoms == 0 )
-	{
-		GEOSSetSRID(geos_result, srid);
-		return geos_result;
-	}
+/* Find holes of each face */
+static void
+findFaceHoles(Face** faces, int nfaces)
+{
+  int i, j, h;
+
+  /* We sort by envelope area so that we know holes are only
+   * after their shells */
+  qsort(faces, nfaces, sizeof(Face*), compare_by_envarea);
+  for (i=0; i<nfaces; ++i) {
+    Face* f = faces[i];
+    int nholes = GEOSGetNumInteriorRings(f->geom);
+    LWDEBUGF(2, "Scanning face %d with env area %g and %d holes", i, f->envarea, nholes);
+    for (h=0; h<nholes; ++h) {
+      const GEOSGeometry *hole = GEOSGetInteriorRingN(f->geom, h);
+      LWDEBUGF(2, "Looking for hole %d/%d of face %d among %d other faces", h+1, nholes, i, nfaces-i-1);
+      for (j=i+1; j<nfaces; ++j) {
+        Face* f2 = faces[j];
+        if ( f2->parent ) continue; /* hole already assigned */
+        /* TODO: short circuit: check envelopes ? */
+        const GEOSGeometry *f2er = GEOSGetExteriorRing(f2->geom); /* TODO: cache this ? */
+        if ( GEOSEquals(f2er, hole) ) {
+          LWDEBUGF(2, "Hole %d/%d of face %d is face %d", h+1, nholes, i, j);
+          f2->parent = f;
+          break;
+        }
+      }
+    }
+  }
+}
 
-	/*
-	 * Return first geometry if we only have one in collection,
-	 * to avoid the unnecessary Geometry clone below.
-	 */
-	if ( ngeoms == 1 )
-	{
-		tmp = (GEOSGeometry *)GEOSGetGeometryN(geos_result, 0);
-		if ( ! tmp )
-		{
-			GEOSGeom_destroy(geos_result);
-			return 0; /* exception */
-		}
-		shp = GEOSGeom_clone(tmp);
-		GEOSGeom_destroy(geos_result); /* only safe after the clone above */
-		GEOSSetSRID(shp, srid);
-		return shp;
-	}
+static GEOSGeometry*
+collectFacesWithEvenParens(Face** faces, int nfaces)
+{
+  GEOSGeometry **geoms = lwalloc(sizeof(GEOSGeometry*)*nfaces);
+  GEOSGeometry *ret;
+  unsigned int ngeoms = 0;
+  int i;
+
+  for (i=0; i<nfaces; ++i) {
+    Face *f = faces[i];
+    if ( countParens(f) % 2 ) continue; /* we skip odd parents geoms */
+    geoms[ngeoms++] = GEOSGeom_clone(f->geom);
+  }
+
+  ret = GEOSGeom_createCollection(GEOS_MULTIPOLYGON, geoms, ngeoms);
+  lwfree(geoms);
+  return ret;
+}
 
-	/*
-	 * Iteratively invoke symdifference on outer rings
-	 * as suggested by Carl Anderson:
-	 * postgis-devel/2005-December/001805.html
-	 *
-	 * Order by number of points, to speed things up.
-	 * See http://trac.osgeo.org/postgis/ticket/1806
-	 */
-	geoms = lwalloc(sizeof(GEOSGeometry*)*ngeoms);
-	for (i=0; i<ngeoms; ++i)
-		geoms[i] = GEOSGetGeometryN(geos_result, i);
-	qsort(geoms, ngeoms, sizeof(GEOSGeometry*), compare_by_numpoints);
-
-	shp = NULL;
-	for (i=0; i<ngeoms; ++i)
-	{
-		GEOSGeom extring;
-		extring = GEOSGeom_createPolygon(
-		              GEOSGeom_clone(GEOSGetExteriorRing(geoms[i])),
-		              NULL, 0
-		          );
+GEOSGeometry*
+LWGEOM_GEOS_buildArea(const GEOSGeometry* geom_in)
+{
+  GEOSGeometry *tmp;
+  GEOSGeometry *geos_result, *shp;
+  GEOSGeometry const *vgeoms[1];
+  uint32_t i, ngeoms;
+  int srid = GEOSGetSRID(geom_in);
+  Face ** geoms;
 
-		if ( extring == NULL ) /* exception */
-		{
-			lwfree(geoms);
-			lwerror("GEOSCreatePolygon threw an exception");
-			return 0;
-		}
+  vgeoms[0] = geom_in;
+  geos_result = GEOSPolygonize(vgeoms, 1);
 
-		if ( shp == NULL )
-		{
-			shp = extring;
-			LWDEBUGF(3, "GEOSpolygonize: shp:%s",
-			               lwgeom_to_ewkt(GEOS2LWGEOM(shp, 0)));
-		}
-		else
-		{
-			tmp = GEOSSymDifference(shp, extring);
-			if ( tmp == NULL ) /* exception */
-			{
-				lwfree(geoms);
-				lwerror("GEOSSymDifference threw an exception: %s", lwgeom_geos_errmsg);
-				return NULL;
-			}
-			LWDEBUGF(3, "GEOSpolygonize: SymDifference(%s, %s):%s",
-			               lwgeom_to_ewkt(GEOS2LWGEOM(shp, 0)),
-			               lwgeom_to_ewkt(GEOS2LWGEOM(extring, 0)),
-			               lwgeom_to_ewkt(GEOS2LWGEOM(tmp, 0))
-			              );
-			GEOSGeom_destroy(shp);
-			GEOSGeom_destroy(extring);
-			shp = tmp;
-		}
-	}
+  LWDEBUGF(3, "GEOSpolygonize returned @ %p", geos_result);
 
-	lwfree(geoms);
-	GEOSGeom_destroy(geos_result);
+  /* Null return from GEOSpolygonize (an exception) */
+  if ( ! geos_result ) return 0;
 
-	GEOSSetSRID(shp, srid);
+  /*
+   * We should now have a collection
+   */
+#if PARANOIA_LEVEL > 0
+  if ( GEOSGeometryTypeId(geos_result) != COLLECTIONTYPE )
+  {
+    GEOSGeom_destroy(geos_result);
+    lwerror("Unexpected return from GEOSpolygonize");
+    return 0;
+  }
+#endif
 
-	return shp;
+  ngeoms = GEOSGetNumGeometries(geos_result);
+
+  LWDEBUGF(3, "GEOSpolygonize: ngeoms in polygonize output: %d", ngeoms);
+  LWDEBUGF(3, "GEOSpolygonize: polygonized:%s",
+              lwgeom_to_ewkt(GEOS2LWGEOM(geos_result, 0)));
+
+  /*
+   * No geometries in collection, early out
+   */
+  if ( ngeoms == 0 )
+  {
+    GEOSSetSRID(geos_result, srid);
+    return geos_result;
+  }
+
+  /*
+   * Return first geometry if we only have one in collection,
+   * to avoid the unnecessary Geometry clone below.
+   */
+  if ( ngeoms == 1 )
+  {
+    tmp = (GEOSGeometry *)GEOSGetGeometryN(geos_result, 0);
+    if ( ! tmp )
+    {
+      GEOSGeom_destroy(geos_result);
+      return 0; /* exception */
+    }
+    shp = GEOSGeom_clone(tmp);
+    GEOSGeom_destroy(geos_result); /* only safe after the clone above */
+    GEOSSetSRID(shp, srid);
+    return shp;
+  }
+
+  LWDEBUGF(2, "Polygonize returned %d geoms", ngeoms);
+
+  /*
+   * Iteratively invoke symdifference on outer rings
+   * as suggested by Carl Anderson:
+   * postgis-devel/2005-December/001805.html
+   *
+   * Order by number of points, to speed things up.
+   * See http://trac.osgeo.org/postgis/ticket/1806
+   */
+  geoms = lwalloc(sizeof(Face**)*ngeoms);
+  for (i=0; i<ngeoms; ++i)
+    geoms[i] = newFace(GEOSGetGeometryN(geos_result, i));
+
+  findFaceHoles(geoms, ngeoms);
+
+  tmp = collectFacesWithEvenParens(geoms, ngeoms);
+
+  for (i=0; i<ngeoms; ++i) delFace(geoms[i]);
+  lwfree(geoms);
+
+  /* geos_result was needed by faces */
+  GEOSGeom_destroy(geos_result);
+
+  shp = GEOSUnionCascaded(tmp);
+  GEOSGeom_destroy(tmp);
+
+
+  GEOSSetSRID(shp, srid);
+
+  return shp;
 }
 
 LWGEOM*
@@ -916,6 +985,8 @@ lwgeom_buildarea(const LWGEOM *geom)
 	return geom_out;
 }
 
+/* ------------ end of BuildArea stuff ---------------------------------------------------------------------} */
+
 LWGEOM*
 lwgeom_geos_noop(const LWGEOM* geom_in)
 {