From b6023cb85cdb6029297af712cdd0dc03237408e6 Mon Sep 17 00:00:00 2001 From: Sandro Santilli Date: Thu, 10 May 2012 07:24:08 +0000 Subject: [PATCH] Cleanup new BuildArea code, add profiling prints (#1806) git-svn-id: http://svn.osgeo.org/postgis/trunk@9732 b70326c6-7e19-0410-871a-916f4a2858ee --- liblwgeom/lwgeom_geos.c | 84 ++++++++++++++++++++++++++++++----- liblwgeom/lwgeom_geos_clean.c | 44 +++++++++++++++++- 2 files changed, 117 insertions(+), 11 deletions(-) diff --git a/liblwgeom/lwgeom_geos.c b/liblwgeom/lwgeom_geos.c index d968f0d4d..5edd3c433 100644 --- a/liblwgeom/lwgeom_geos.c +++ b/liblwgeom/lwgeom_geos.c @@ -17,6 +17,8 @@ #include +#undef LWGEOM_PROFILE_BUILDAREA + #define LWGEOM_GEOS_ERRMSG_MAXSIZE 256 char lwgeom_geos_errmsg[LWGEOM_GEOS_ERRMSG_MAXSIZE]; @@ -799,8 +801,12 @@ findFaceHoles(Face** faces, int nfaces) for (j=i+1; jparent ) continue; /* hole already assigned */ - /* TODO: short circuit: check envelopes ? */ - const GEOSGeometry *f2er = GEOSGetExteriorRing(f2->geom); /* TODO: cache this ? */ + const GEOSGeometry *f2er = GEOSGetExteriorRing(f2->geom); + /* TODO: can be optimized as the ring would have the + * same vertices, possibly in different order. + * maybe comparing number of points could already be + * useful. + */ if ( GEOSEquals(f2er, hole) ) { LWDEBUGF(2, "Hole %d/%d of face %d is face %d", h+1, nholes, i, j); f2->parent = f; @@ -812,7 +818,7 @@ findFaceHoles(Face** faces, int nfaces) } static GEOSGeometry* -collectFacesWithEvenParens(Face** faces, int nfaces) +collectFacesWithEvenAncestors(Face** faces, int nfaces) { GEOSGeometry **geoms = lwalloc(sizeof(GEOSGeometry*)*nfaces); GEOSGeometry *ret; @@ -841,6 +847,9 @@ LWGEOM_GEOS_buildArea(const GEOSGeometry* geom_in) Face ** geoms; vgeoms[0] = geom_in; +#ifdef LWGEOM_PROFILE_BUILDAREA + lwnotice("Polygonizing"); +#endif geos_result = GEOSPolygonize(vgeoms, 1); LWDEBUGF(3, "GEOSpolygonize returned @ %p", geos_result); @@ -861,6 +870,10 @@ LWGEOM_GEOS_buildArea(const GEOSGeometry* geom_in) #endif ngeoms = GEOSGetNumGeometries(geos_result); +#ifdef LWGEOM_PROFILE_BUILDAREA + lwnotice("Num geometries from polygonizer: %d", ngeoms); +#endif + LWDEBUGF(3, "GEOSpolygonize: ngeoms in polygonize output: %d", ngeoms); LWDEBUGF(3, "GEOSpolygonize: polygonized:%s", @@ -896,30 +909,81 @@ LWGEOM_GEOS_buildArea(const GEOSGeometry* geom_in) LWDEBUGF(2, "Polygonize returned %d geoms", ngeoms); /* - * Iteratively invoke symdifference on outer rings - * as suggested by Carl Anderson: - * postgis-devel/2005-December/001805.html + * Polygonizer returns a polygon for each face in the built topology. + * + * This means that for any face with holes we'll have other faces + * representing each hole. We can imagine a parent-child relationship + * between these faces. + * + * In order to maximize the number of visible rings in output we + * only use those faces which have an even number of parents. * - * Order by number of points, to speed things up. + * Example: + * + * +---------------+ + * | L0 | L0 has no parents + * | +---------+ | + * | | L1 | | L1 is an hole of L0 + * | | +---+ | | + * | | |L2 | | | L2 is an hole of L1 (which is an hole of L0) + * | | | | | | + * | | +---+ | | + * | +---------+ | + * | | + * +---------------+ + * * See http://trac.osgeo.org/postgis/ticket/1806 + * */ + +#ifdef LWGEOM_PROFILE_BUILDAREA + lwnotice("Preparing face structures"); +#endif + + /* Prepare face structures for later analysis */ geoms = lwalloc(sizeof(Face**)*ngeoms); for (i=0; i /* #define POSTGIS_DEBUG_LEVEL 4 */ +#undef LWGEOM_PROFILE_MAKEVALID /* @@ -409,6 +410,10 @@ LWGEOM_GEOS_makeValidPolygon(const GEOSGeometry* gin) /* Use noded boundaries as initial "cut" edges */ +#ifdef LWGEOM_PROFILE_MAKEVALID + lwnotice("ST_MakeValid: noding lines"); +#endif + geos_cut_edges = LWGEOM_GEOS_nodeLines(geos_bound); if ( NULL == geos_cut_edges ) @@ -424,6 +429,10 @@ LWGEOM_GEOS_makeValidPolygon(const GEOSGeometry* gin) GEOSGeometry* pi; GEOSGeometry* po; +#ifdef LWGEOM_PROFILE_MAKEVALID + lwnotice("ST_MakeValid: extracting unique points from bounds"); +#endif + pi = GEOSGeom_extractUniquePoints(geos_bound); if ( NULL == pi ) { @@ -437,6 +446,10 @@ LWGEOM_GEOS_makeValidPolygon(const GEOSGeometry* gin) "Boundaries input points %s", lwgeom_to_ewkt(GEOS2LWGEOM(pi, 0))); +#ifdef LWGEOM_PROFILE_MAKEVALID + lwnotice("ST_MakeValid: extracting unique points from cut_edges"); +#endif + po = GEOSGeom_extractUniquePoints(geos_cut_edges); if ( NULL == po ) { @@ -451,6 +464,10 @@ LWGEOM_GEOS_makeValidPolygon(const GEOSGeometry* gin) "Boundaries output points %s", lwgeom_to_ewkt(GEOS2LWGEOM(po, 0))); +#ifdef LWGEOM_PROFILE_MAKEVALID + lwnotice("ST_MakeValid: find collapse points"); +#endif + collapse_points = GEOSDifference(pi, po); if ( NULL == collapse_points ) { @@ -465,6 +482,10 @@ LWGEOM_GEOS_makeValidPolygon(const GEOSGeometry* gin) "Collapse points: %s", lwgeom_to_ewkt(GEOS2LWGEOM(collapse_points, 0))); +#ifdef LWGEOM_PROFILE_MAKEVALID + lwnotice("ST_MakeValid: cleanup(1)"); +#endif + GEOSGeom_destroy(pi); GEOSGeom_destroy(po); } @@ -483,7 +504,6 @@ LWGEOM_GEOS_makeValidPolygon(const GEOSGeometry* gin) return NULL; } - /* * See if an area can be build with the remaining edges * and if it can, symdifference with the original area. @@ -497,6 +517,10 @@ LWGEOM_GEOS_makeValidPolygon(const GEOSGeometry* gin) GEOSGeometry* symdif=0; GEOSGeometry* new_cut_edges=0; +#ifdef LWGEOM_PROFILE_MAKEVALID + lwnotice("ST_MakeValid: building area from %d edges", GEOSGetNumGeometries(geos_cut_edges)); +#endif + /* * ASSUMPTION: cut_edges should already be fully noded */ @@ -522,6 +546,9 @@ LWGEOM_GEOS_makeValidPolygon(const GEOSGeometry* gin) * We succeeded in building a ring ! */ +#ifdef LWGEOM_PROFILE_MAKEVALID + lwnotice("ST_MakeValid: ring built with %d cut edges, saving boundaries", GEOSGetNumGeometries(geos_cut_edges)); +#endif /* * Save the new ring boundaries first (to compute @@ -540,6 +567,10 @@ LWGEOM_GEOS_makeValidPolygon(const GEOSGeometry* gin) return NULL; } +#ifdef LWGEOM_PROFILE_MAKEVALID + lwnotice("ST_MakeValid: running SymDifference with new area"); +#endif + /* * Now symdif new and old area */ @@ -566,8 +597,15 @@ LWGEOM_GEOS_makeValidPolygon(const GEOSGeometry* gin) * ASSUMPTION: only the previous cut-edges can be * left, so we don't need to reconsider * the whole original boundaries + * + * NOTE: this is an expensive operation. + * */ +#ifdef LWGEOM_PROFILE_MAKEVALID + lwnotice("ST_MakeValid: computing new cut_edges (GEOSDifference)"); +#endif + new_cut_edges = GEOSDifference(geos_cut_edges, new_area_bound); GEOSGeom_destroy(new_area_bound); if ( ! new_cut_edges ) /* an exception ? */ @@ -583,6 +621,10 @@ LWGEOM_GEOS_makeValidPolygon(const GEOSGeometry* gin) geos_cut_edges = new_cut_edges; } +#ifdef LWGEOM_PROFILE_MAKEVALID + lwnotice("ST_MakeValid: final checks"); +#endif + if ( ! GEOSisEmpty(geos_area) ) { vgeoms[nvgeoms++] = geos_area; -- 2.40.0