From 75a2406a9f3a4812b28f524e083bd80951c33b5e Mon Sep 17 00:00:00 2001 From: Sandro Santilli Date: Sat, 30 Sep 2017 15:26:28 +0000 Subject: [PATCH] Add lwt_Polygonize and lwt_AddLineNoFace functions This were initially written for librttopo. The plan in PostGIS is to use them for ST_CreateTopoGeo, which already needs to load all input geometry in memory. git-svn-id: http://svn.osgeo.org/postgis/trunk@15857 b70326c6-7e19-0410-871a-916f4a2858ee --- liblwgeom/liblwgeom_topo.h | 56 +- liblwgeom/lwgeom_topo.c | 1311 +++++++++++++++++++++++++++++++----- 2 files changed, 1199 insertions(+), 168 deletions(-) diff --git a/liblwgeom/liblwgeom_topo.h b/liblwgeom/liblwgeom_topo.h index 183af5b0f..ceb2cda3f 100644 --- a/liblwgeom/liblwgeom_topo.h +++ b/liblwgeom/liblwgeom_topo.h @@ -18,7 +18,7 @@ * ********************************************************************** * - * Copyright (C) 2015 Sandro Santilli + * Copyright (C) 2015-2017 Sandro Santilli * **********************************************************************/ @@ -1056,6 +1056,60 @@ LWT_ELEMID lwt_AddPoint(LWT_TOPOLOGY* topo, LWPOINT* point, double tol); LWT_ELEMID* lwt_AddLine(LWT_TOPOLOGY* topo, LWLINE* line, double tol, int* nedges); +/** + * Adds a linestring to the topology without determining generated faces + * + * The given line will snap to existing nodes or edges within given + * tolerance. Existing edges or faces may be split by the line. + * + * Side faces for the new edges will not be determined and no new + * faces will be created, effectively leaving the topology in an + * invalid state (WARNING!) + * + * @param topo the topology to operate on + * @param line the line to add + * @param tol snap tolerance, the topology tolerance will be used if 0 + * @param nedges output parameter, will be set to number of edges the + * line was split into, or -1 on error + * (liblwgeom error handler will be invoked with error message) + * + * @return an array of edge identifiers that sewed togheter + * will build up the input linestring (after snapping). Caller + * will need to free the array using lwfree(), if not null. + */ +LWT_ELEMID* lwt_AddLineNoFace(LWT_TOPOLOGY* topo, LWLINE* line, double tol, + int* nedges); + +/* + * Determine and register all topology faces: + * + * - Determines which faces are generated by existing + * edges. + * - Creates face records with correct mbr + * - Update edge left/right face attributes + * + * Precondition: + * - the topology edges are correctly linked + * - there are no faces registered in the topology + * + * Postconditions: + * - all left/right face attributes of edges + * reference faces with correct mbr. + * + * Notes: + * - does not attempt to assign isolated nodes to their + * containing faces + * - does not remove existing face records + * - loads in memory all the topology edges + * + * @param topo the topology to operate on + * + * @return 0 on success, -1 on error + * (librtgeom error handler will be invoked with error + * message) + */ +int lwt_Polygonize(LWT_TOPOLOGY* topo); + /** * Adds a polygon to the topology * diff --git a/liblwgeom/lwgeom_topo.c b/liblwgeom/lwgeom_topo.c index 17535cfac..ba48c60b9 100644 --- a/liblwgeom/lwgeom_topo.c +++ b/liblwgeom/lwgeom_topo.c @@ -26,7 +26,7 @@ #include "../postgis_config.h" -/*#define POSTGIS_DEBUG_LEVEL 1*/ +#define POSTGIS_DEBUG_LEVEL 1 #include "lwgeom_log.h" #include "liblwgeom_internal.h" @@ -521,9 +521,17 @@ lwt_FreeTopology( LWT_TOPOLOGY* topo ) lwfree(topo); } -LWT_ELEMID -lwt_AddIsoNode( LWT_TOPOLOGY* topo, LWT_ELEMID face, - LWPOINT* pt, int skipISOChecks ) +/** + * @param checkFace if non zero will check the given face + * for really containing the point or determine the + * face when given face is -1. Use 0 to simply use + * the given face value, no matter what (effectively + * allowing adding a non-isolated point when used + * with face=-1). + */ +static LWT_ELEMID +_lwt_AddIsoNode( LWT_TOPOLOGY* topo, LWT_ELEMID face, + LWPOINT* pt, int skipISOChecks, int checkFace ) { LWT_ELEMID foundInFace = -1; @@ -541,7 +549,7 @@ lwt_AddIsoNode( LWT_TOPOLOGY* topo, LWT_ELEMID face, } } - if ( face == -1 || ! skipISOChecks ) + if ( checkFace && ( face == -1 || ! skipISOChecks ) ) { foundInFace = lwt_be_getFaceContainingPoint(topo, pt); /*x*/ if ( foundInFace == -2 ) { @@ -577,6 +585,13 @@ lwt_AddIsoNode( LWT_TOPOLOGY* topo, LWT_ELEMID face, return node.node_id; } +LWT_ELEMID +lwt_AddIsoNode( LWT_TOPOLOGY* topo, LWT_ELEMID face, + LWPOINT* pt, int skipISOChecks ) +{ + return _lwt_AddIsoNode( topo, face, pt, skipISOChecks, 1 ); +} + /* Check that an edge does not cross an existing node or edge * * @param myself the id of an edge to skip, if any @@ -1754,61 +1769,12 @@ _lwt_GetInteriorEdgePoint(const LWLINE* edge, POINT2D* ip) return 1; } -/* - * Add a split face by walking on the edge side. - * - * @param topo the topology to act upon - * @param sedge edge id and walking side and direction - * (forward,left:positive backward,right:negative) - * @param face the face in which the edge identifier is known to be - * @param mbr_only do not create a new face but update MBR of the current - * - * @return: - * -1: if mbr_only was requested - * 0: if the edge does not form a ring - * -1: if it is impossible to create a face on the requested side - * ( new face on the side is the universe ) - * -2: error - * >0 : id of newly added face - */ -static LWT_ELEMID -_lwt_AddFaceSplit( LWT_TOPOLOGY* topo, - LWT_ELEMID sedge, LWT_ELEMID face, - int mbr_only ) +static LWPOLY * +_lwt_MakeRingShell(LWT_TOPOLOGY *topo, LWT_ELEMID *signed_edge_ids, int num_signed_edge_ids) { - int numedges, numfaceedges, i, j; - int newface_outside; - int num_signed_edge_ids; - LWT_ELEMID *signed_edge_ids; LWT_ELEMID *edge_ids; - LWT_ISO_EDGE *edges; + int numedges, i, j; LWT_ISO_EDGE *ring_edges; - LWT_ISO_EDGE *forward_edges = NULL; - int forward_edges_count = 0; - LWT_ISO_EDGE *backward_edges = NULL; - int backward_edges_count = 0; - - signed_edge_ids = lwt_be_getRingEdges(topo, sedge, - &num_signed_edge_ids, 0); - if ( ! signed_edge_ids ) { - lwerror("Backend error (no ring edges for edge %" LWTFMT_ELEMID "): %s", - sedge, lwt_be_lastErrorMessage(topo->be_iface)); - return -2; - } - LWDEBUGF(1, "getRingEdges returned %d edges", num_signed_edge_ids); - - /* You can't get to the other side of an edge forming a ring */ - for (i=0; ibe_iface)); - return -2; + return NULL; } else if ( i != numedges ) { lwfree( signed_edge_ids ); _lwt_release_edges(ring_edges, numedges); lwerror("Unexpected error: %d edges found when expecting %d", i, numedges); - return -2; + return NULL; } /* Should now build a polygon with those edges, in the order @@ -1851,8 +1815,7 @@ _lwt_AddFaceSplit( LWT_TOPOLOGY* topo, for ( i=0; i0 : id of newly added face + */ +static LWT_ELEMID +_lwt_AddFaceSplit( LWT_TOPOLOGY* topo, + LWT_ELEMID sedge, LWT_ELEMID face, + int mbr_only ) +{ + int numfaceedges, i, j; + int newface_outside; + int num_signed_edge_ids; + LWT_ELEMID *signed_edge_ids; + LWT_ISO_EDGE *edges; + LWT_ISO_EDGE *forward_edges = NULL; + int forward_edges_count = 0; + LWT_ISO_EDGE *backward_edges = NULL; + int backward_edges_count = 0; + + signed_edge_ids = lwt_be_getRingEdges(topo, sedge, + &num_signed_edge_ids, 0); + if ( ! signed_edge_ids ) { + lwerror("Backend error (no ring edges for edge %" LWTFMT_ELEMID "): %s", + sedge, lwt_be_lastErrorMessage(topo->be_iface)); + return -2; + } + LWDEBUGF(1, "getRingEdges returned %d edges", num_signed_edge_ids); + + /* You can't get to the other side of an edge forming a ring */ + for (i=0; ibe_iface)); + return -2; + } + const POINTARRAY *pa = shell->rings[0]; int isccw = ptarray_isccw(pa); LWDEBUGF(1, "Ring of edge %" LWTFMT_ELEMID " is %sclockwise", @@ -1911,7 +1942,6 @@ _lwt_AddFaceSplit( LWT_TOPOLOGY* topo, { lwpoly_free(shell); lwfree( signed_edge_ids ); - _lwt_release_edges(ring_edges, numedges); /* Face on the left side of this ring is the universe face. * Next call (for the other side) should create the split face */ @@ -1932,7 +1962,6 @@ _lwt_AddFaceSplit( LWT_TOPOLOGY* topo, if ( ret == -1 ) { lwfree( signed_edge_ids ); - _lwt_release_edges(ring_edges, numedges); lwpoly_free(shell); /* NOTE: owns shellbox above */ lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -2; @@ -1940,14 +1969,12 @@ _lwt_AddFaceSplit( LWT_TOPOLOGY* topo, if ( ret != 1 ) { lwfree( signed_edge_ids ); - _lwt_release_edges(ring_edges, numedges); lwpoly_free(shell); /* NOTE: owns shellbox above */ lwerror("Unexpected error: %d faces found when expecting 1", ret); return -2; } }} lwfree( signed_edge_ids ); - _lwt_release_edges(ring_edges, numedges); lwpoly_free(shell); /* NOTE: owns shellbox above */ return -1; /* mbr only was requested */ } @@ -1964,7 +1991,6 @@ _lwt_AddFaceSplit( LWT_TOPOLOGY* topo, { lwfree( signed_edge_ids ); lwpoly_free(shell); /* NOTE: owns shellbox */ - _lwt_release_edges(ring_edges, numedges); lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -2; } @@ -1972,7 +1998,6 @@ _lwt_AddFaceSplit( LWT_TOPOLOGY* topo, { lwfree( signed_edge_ids ); lwpoly_free(shell); /* NOTE: owns shellbox */ - _lwt_release_edges(ring_edges, numedges); lwerror("Unexpected error: %d faces found when expecting 1", nfaces); return -2; } @@ -1989,7 +2014,6 @@ _lwt_AddFaceSplit( LWT_TOPOLOGY* topo, { lwfree( signed_edge_ids ); lwpoly_free(shell); /* NOTE: owns shellbox */ - _lwt_release_edges(ring_edges, numedges); lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -2; } @@ -1997,7 +2021,6 @@ _lwt_AddFaceSplit( LWT_TOPOLOGY* topo, { lwfree( signed_edge_ids ); lwpoly_free(shell); /* NOTE: owns shellbox */ - _lwt_release_edges(ring_edges, numedges); lwerror("Unexpected error: %d faces inserted when expecting 1", ret); return -2; } @@ -2031,7 +2054,6 @@ _lwt_AddFaceSplit( LWT_TOPOLOGY* topo, edges = lwt_be_getEdgeByFace( topo, &face, &numfaceedges, fields, newface.mbr ); if ( numfaceedges == -1 ) { lwfree( signed_edge_ids ); - _lwt_release_edges(ring_edges, numedges); lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -2; } @@ -2042,7 +2064,6 @@ _lwt_AddFaceSplit( LWT_TOPOLOGY* topo, if ( ! shellgg ) { lwpoly_free(shell); lwfree(signed_edge_ids); - _lwt_release_edges(ring_edges, numedges); _lwt_release_edges(edges, numfaceedges); lwerror("Could not convert shell geometry to GEOS: %s", lwgeom_geos_errmsg); return -2; @@ -2052,7 +2073,6 @@ _lwt_AddFaceSplit( LWT_TOPOLOGY* topo, GEOSGeom_destroy(shellgg); lwpoly_free(shell); lwfree(signed_edge_ids); - _lwt_release_edges(ring_edges, numedges); _lwt_release_edges(edges, numfaceedges); lwerror("Could not prepare shell geometry: %s", lwgeom_geos_errmsg); return -2; @@ -2111,9 +2131,8 @@ _lwt_AddFaceSplit( LWT_TOPOLOGY* topo, GEOSGeom_destroy(shellgg); lwfree(signed_edge_ids); lwpoly_free(shell); - lwfree(forward_edges); /* contents owned by ring_edges */ - lwfree(backward_edges); /* contents owned by ring_edges */ - _lwt_release_edges(ring_edges, numedges); + lwfree(forward_edges); /* contents owned by edges */ + lwfree(backward_edges); /* contents owned by edges */ _lwt_release_edges(edges, numfaceedges); lwerror("Could not find interior point for edge %d: %s", e->edge_id, lwgeom_geos_errmsg); @@ -2128,9 +2147,8 @@ _lwt_AddFaceSplit( LWT_TOPOLOGY* topo, GEOSGeom_destroy(shellgg); lwfree(signed_edge_ids); lwpoly_free(shell); - lwfree(forward_edges); /* contents owned by ring_edges */ - lwfree(backward_edges); /* contents owned by ring_edges */ - _lwt_release_edges(ring_edges, numedges); + lwfree(forward_edges); /* contents owned by edges */ + lwfree(backward_edges); /* contents owned by edges */ _lwt_release_edges(edges, numfaceedges); lwerror("Could not convert edge geometry to GEOS: %s", lwgeom_geos_errmsg); @@ -2148,9 +2166,8 @@ _lwt_AddFaceSplit( LWT_TOPOLOGY* topo, GEOSGeom_destroy(shellgg); lwfree(signed_edge_ids); lwpoly_free(shell); - lwfree(forward_edges); /* contents owned by ring_edges */ - lwfree(backward_edges); /* contents owned by ring_edges */ - _lwt_release_edges(ring_edges, numedges); + lwfree(forward_edges); /* contents owned by edges */ + lwfree(backward_edges); /* contents owned by edges */ _lwt_release_edges(edges, numfaceedges); lwerror("GEOS exception on PreparedContains: %s", lwgeom_geos_errmsg); return -2; @@ -2242,7 +2259,6 @@ _lwt_AddFaceSplit( LWT_TOPOLOGY* topo, } - _lwt_release_edges(ring_edges, numedges); _lwt_release_edges(edges, numfaceedges); /* Update isolated nodes which are now in new face */ @@ -2335,6 +2351,13 @@ _lwt_AddFaceSplit( LWT_TOPOLOGY* topo, return newface.face_id; } +/** + * @param modFace can be + * 0 - have two new faces replace a splitted face + * 1 - modify a splitted face, adding a new one + * -1 - do not check at all for face splitting + * + */ static LWT_ELEMID _lwt_AddEdge( LWT_TOPOLOGY* topo, LWT_ELEMID start_node, LWT_ELEMID end_node, @@ -2565,7 +2588,7 @@ _lwt_AddEdge( LWT_TOPOLOGY* topo, newedge.edge_id, newedge.next_left, prev_right); if ( newedge.face_right == -1 ) { newedge.face_right = span.ccwFace; - } else if ( newedge.face_right != epan.ccwFace ) { + } else if ( modFace != -1 && newedge.face_right != epan.ccwFace ) { /* side-location conflict */ lwerror("Side-location conflict: " "new edge starts in face" @@ -2577,7 +2600,7 @@ _lwt_AddEdge( LWT_TOPOLOGY* topo, } if ( newedge.face_left == -1 ) { newedge.face_left = span.cwFace; - } else if ( newedge.face_left != epan.cwFace ) { + } else if ( modFace != -1 && newedge.face_left != epan.cwFace ) { /* side-location conflict */ lwerror("Side-location conflict: " "new edge starts in face" @@ -2608,7 +2631,7 @@ _lwt_AddEdge( LWT_TOPOLOGY* topo, newedge.face_left, newedge.face_right); return -1; } - else if ( newedge.face_left == -1 ) + else if ( newedge.face_left == -1 && modFace > -1 ) { lwerror("Could not derive edge face from linked primitives:" " invalid topology ?"); @@ -2720,7 +2743,9 @@ _lwt_AddEdge( LWT_TOPOLOGY* topo, } } - /* Check face splitting */ + /* Check face splitting, if required */ + + if ( modFace > -1 ) { if ( ! isclosed && ( epan.was_isolated || span.was_isolated ) ) { @@ -2788,6 +2813,8 @@ _lwt_AddEdge( LWT_TOPOLOGY* topo, } } + } // end of face split checking + return newedge.edge_id; } @@ -3601,75 +3628,82 @@ lwt_ChangeEdgeGeom(LWT_TOPOLOGY* topo, LWT_ELEMID edge_id, LWLINE *geom) -- would be larger than the old face MBR... -- */ - int facestoupdate = 0; - LWT_ISO_FACE faces[2]; - LWGEOM *nface1 = NULL; - LWGEOM *nface2 = NULL; - if ( oldedge->face_left != 0 ) - { - nface1 = lwt_GetFaceGeometry(topo, oldedge->face_left); - if ( ! nface1 ) - { - lwerror("lwt_ChangeEdgeGeom could not construct face %" - LWTFMT_ELEMID ", on the left of edge %" LWTFMT_ELEMID, - oldedge->face_left, edge_id); - return -1; - } -#if 0 - { - size_t sz; - char *wkt = lwgeom_to_wkt(nface1, WKT_EXTENDED, 2, &sz); - LWDEBUGF(1, "new geometry of face left (%d): %s", (int)oldedge->face_left, wkt); - lwfree(wkt); - } -#endif - lwgeom_add_bbox(nface1); - faces[facestoupdate].face_id = oldedge->face_left; - /* ownership left to nface */ - faces[facestoupdate++].mbr = nface1->bbox; - } - if ( oldedge->face_right != 0 - /* no need to update twice the same face.. */ - && oldedge->face_right != oldedge->face_left ) - { - nface2 = lwt_GetFaceGeometry(topo, oldedge->face_right); - if ( ! nface2 ) + const GBOX* oldbox = lwgeom_get_bbox(lwline_as_lwgeom(oldedge->geom)); + const GBOX* newbox = lwgeom_get_bbox(lwline_as_lwgeom(geom)); + if ( ! gbox_same(oldbox, newbox) ) + {{ + int facestoupdate = 0; + LWT_ISO_FACE faces[2]; + LWGEOM *nface1 = NULL; + LWGEOM *nface2 = NULL; + if ( oldedge->face_left > 0 ) { - lwerror("lwt_ChangeEdgeGeom could not construct face %" - LWTFMT_ELEMID ", on the right of edge %" LWTFMT_ELEMID, - oldedge->face_right, edge_id); - return -1; - } -#if 0 + nface1 = lwt_GetFaceGeometry(topo, oldedge->face_left); + if ( ! nface1 ) + { + lwerror("lwt_ChangeEdgeGeom could not construct face %" + LWTFMT_ELEMID ", on the left of edge %" LWTFMT_ELEMID, + oldedge->face_left, edge_id); + return -1; + } + #if 0 + { + size_t sz; + char *wkt = lwgeom_to_wkt(nface1, WKT_EXTENDED, 2, &sz); + LWDEBUGF(1, "new geometry of face left (%d): %s", (int)oldedge->face_left, wkt); + lwfree(wkt); + } + #endif + lwgeom_add_bbox(nface1); + faces[facestoupdate].face_id = oldedge->face_left; + /* ownership left to nface */ + faces[facestoupdate++].mbr = nface1->bbox; + } + if ( oldedge->face_right > 0 + /* no need to update twice the same face.. */ + && oldedge->face_right != oldedge->face_left ) { - size_t sz; - char *wkt = lwgeom_to_wkt(nface2, WKT_EXTENDED, 2, &sz); - LWDEBUGF(1, "new geometry of face right (%d): %s", (int)oldedge->face_right, wkt); - lwfree(wkt); + nface2 = lwt_GetFaceGeometry(topo, oldedge->face_right); + if ( ! nface2 ) + { + lwerror("lwt_ChangeEdgeGeom could not construct face %" + LWTFMT_ELEMID ", on the right of edge %" LWTFMT_ELEMID, + oldedge->face_right, edge_id); + return -1; + } + #if 0 + { + size_t sz; + char *wkt = lwgeom_to_wkt(nface2, WKT_EXTENDED, 2, &sz); + LWDEBUGF(1, "new geometry of face right (%d): %s", (int)oldedge->face_right, wkt); + lwfree(wkt); + } + #endif + lwgeom_add_bbox(nface2); + faces[facestoupdate].face_id = oldedge->face_right; + faces[facestoupdate++].mbr = nface2->bbox; /* ownership left to nface */ } -#endif - lwgeom_add_bbox(nface2); - faces[facestoupdate].face_id = oldedge->face_right; - faces[facestoupdate++].mbr = nface2->bbox; /* ownership left to nface */ - } - LWDEBUGF(1, "%d faces to update", facestoupdate); - if ( facestoupdate ) - { - i = lwt_be_updateFacesById( topo, &(faces[0]), facestoupdate ); - if ( i != facestoupdate ) + LWDEBUGF(1, "%d faces to update", facestoupdate); + if ( facestoupdate ) { - if ( nface1 ) lwgeom_free(nface1); - if ( nface2 ) lwgeom_free(nface2); - _lwt_release_edges(oldedge, 1); - if ( i == -1 ) - lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); - else - lwerror("Unexpected error: %d faces found when expecting 1", i); - return -1; + i = lwt_be_updateFacesById( topo, &(faces[0]), facestoupdate ); + if ( i != facestoupdate ) + { + if ( nface1 ) lwgeom_free(nface1); + if ( nface2 ) lwgeom_free(nface2); + _lwt_release_edges(oldedge, 1); + if ( i == -1 ) + lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); + else + lwerror("Unexpected error: %d faces found when expecting 1", i); + return -1; + } } - } - if ( nface1 ) lwgeom_free(nface1); - if ( nface2 ) lwgeom_free(nface2); + if ( nface1 ) lwgeom_free(nface1); + if ( nface2 ) lwgeom_free(nface2); + }} else {{ + lwnotice("BBOX of edge changed edge did not change"); + }} LWDEBUG(1, "all done, cleaning up edges"); @@ -5017,8 +5051,13 @@ compare_scored_pointer(const void *si1, const void *si2) return 0; } -LWT_ELEMID -lwt_AddPoint(LWT_TOPOLOGY* topo, LWPOINT* point, double tol) +/* + * @param findFace if non-zero the code will determine which face + * contains the given point (unless it is known to be NOT + * isolated) + */ +static LWT_ELEMID +_lwt_AddPoint(LWT_TOPOLOGY* topo, LWPOINT* point, double tol, int findFace) { int num, i; double mindist = FLT_MAX; @@ -5077,7 +5116,8 @@ lwt_AddPoint(LWT_TOPOLOGY* topo, LWPOINT* point, double tol) LWGEOM *g = lwpoint_as_lwgeom(n->geom); double dist = lwgeom_mindistance2d(g, pt); /* TODO: move this check in the previous sort scan ... */ - if ( dist >= tol ) continue; /* must be closer than tolerated */ + /* must be closer than tolerated, unless distance is zero */ + if ( dist && dist >= tol ) continue; if ( ! id || dist < mindist ) { id = n->node_id; @@ -5338,7 +5378,7 @@ lwt_AddPoint(LWT_TOPOLOGY* topo, LWPOINT* point, double tol) { /* The point is isolated, add it as such */ /* TODO: pass 1 as last argument (skipChecks) ? */ - id = lwt_AddIsoNode(topo, -1, point, 0); + id = _lwt_AddIsoNode(topo, -1, point, 0, findFace); if ( -1 == id ) { /* should have invoked lwerror already, leaking memory */ @@ -5350,6 +5390,12 @@ lwt_AddPoint(LWT_TOPOLOGY* topo, LWPOINT* point, double tol) return id; } +LWT_ELEMID +lwt_AddPoint(LWT_TOPOLOGY* topo, LWPOINT* point, double tol) +{ + return _lwt_AddPoint(topo, point, tol, 1); +} + /* Return identifier of an equal edge, 0 if none or -1 on error * (and lwerror gets called on error) */ @@ -5421,9 +5467,15 @@ _lwt_GetEqualEdge( LWT_TOPOLOGY *topo, LWLINE *edge ) /* * Add a pre-noded pre-split line edge. Used by lwt_AddLine * Return edge id, 0 if none added (empty edge), -1 on error + * + * @param handleFaceSplit if non-zero the code will check + * if the newly added edge would split a face and if so + * would create new faces accordingly. Otherwise it will + * set left_face and right_face to null (-1) */ static LWT_ELEMID -_lwt_AddLineEdge( LWT_TOPOLOGY* topo, LWLINE* edge, double tol ) +_lwt_AddLineEdge( LWT_TOPOLOGY* topo, LWLINE* edge, double tol, + int handleFaceSplit ) { LWCOLLECTION *col; LWPOINT *start_point, *end_point; @@ -5443,7 +5495,7 @@ _lwt_AddLineEdge( LWT_TOPOLOGY* topo, LWLINE* edge, double tol ) lwnotice("Empty component of noded line"); return 0; /* must be empty */ } - nid[0] = lwt_AddPoint( topo, start_point, tol ); + nid[0] = _lwt_AddPoint( topo, start_point, tol, handleFaceSplit ); lwpoint_free(start_point); /* too late if lwt_AddPoint calls lwerror */ if ( nid[0] == -1 ) return -1; /* lwerror should have been called */ @@ -5454,7 +5506,7 @@ _lwt_AddLineEdge( LWT_TOPOLOGY* topo, LWLINE* edge, double tol ) "after successfully getting first point !?"); return -1; } - nid[1] = lwt_AddPoint( topo, end_point, tol ); + nid[1] = _lwt_AddPoint( topo, end_point, tol, handleFaceSplit ); lwpoint_free(end_point); /* too late if lwt_AddPoint calls lwerror */ if ( nid[1] == -1 ) return -1; /* lwerror should have been called */ @@ -5580,7 +5632,7 @@ _lwt_AddLineEdge( LWT_TOPOLOGY* topo, LWLINE* edge, double tol ) /* TODO: skip checks ? */ - id = lwt_AddEdgeModFace( topo, nid[0], nid[1], edge, 0 ); + id = _lwt_AddEdge( topo, nid[0], nid[1], edge, 0, handleFaceSplit ? 1 : -1 ); LWDEBUGF(1, "lwt_AddEdgeModFace returned %" LWTFMT_ELEMID, id); if ( id == -1 ) { @@ -5616,8 +5668,9 @@ _lwt_split_by_nodes(const LWGEOM *g, const LWGEOM *nodes) return bg; } -LWT_ELEMID* -lwt_AddLine(LWT_TOPOLOGY* topo, LWLINE* line, double tol, int* nedges) +static LWT_ELEMID* +_lwt_AddLine(LWT_TOPOLOGY* topo, LWLINE* line, double tol, int* nedges, + int handleFaceSplit) { LWGEOM *geomsbuf[1]; LWGEOM **geoms; @@ -5680,7 +5733,8 @@ lwt_AddLine(LWT_TOPOLOGY* topo, LWLINE* line, double tol, int* nedges) LWGEOM *g = lwline_as_lwgeom(e->geom); LWDEBUGF(2, "Computing distance from edge %d having %d points", i, e->geom->points->npoints); double dist = lwgeom_mindistance2d(g, noded); - if ( dist >= tol ) continue; /* must be closer than tolerated */ + /* must be closer than tolerated, unless distance is zero */ + if ( dist && dist >= tol ) continue; nearby[nn++] = g; } LWDEBUGF(2, "Found %d lines closer than tolerance (%g)", nn, tol); @@ -5754,7 +5808,8 @@ lwt_AddLine(LWT_TOPOLOGY* topo, LWLINE* line, double tol, int* nedges) LWT_ISO_NODE *n = &(nodes[i]); LWGEOM *g = lwpoint_as_lwgeom(n->geom); double dist = lwgeom_mindistance2d(g, noded); - if ( dist >= tol ) continue; /* must be closer than tolerated */ + /* must be closer than tolerated, unless distance is zero */ + if ( dist && dist >= tol ) continue; nearby[nn++] = g; } if ( nn ) @@ -5843,7 +5898,7 @@ lwt_AddLine(LWT_TOPOLOGY* topo, LWLINE* line, double tol, int* nedges) } #endif - id = _lwt_AddLineEdge( topo, lwgeom_as_lwline(g), tol ); + id = _lwt_AddLineEdge( topo, lwgeom_as_lwline(g), tol, handleFaceSplit ); LWDEBUGF(1, "_lwt_AddLineEdge returned %" LWTFMT_ELEMID, id); if ( id < 0 ) { @@ -5871,6 +5926,18 @@ lwt_AddLine(LWT_TOPOLOGY* topo, LWLINE* line, double tol, int* nedges) return ids; } +LWT_ELEMID* +lwt_AddLine(LWT_TOPOLOGY* topo, LWLINE* line, double tol, int* nedges) +{ + return _lwt_AddLine(topo, line, tol, nedges, 1); +} + +LWT_ELEMID* +lwt_AddLineNoFace(LWT_TOPOLOGY* topo, LWLINE* line, double tol, int* nedges) +{ + return _lwt_AddLine(topo, line, tol, nedges, 0); +} + LWT_ELEMID* lwt_AddPolygon(LWT_TOPOLOGY* topo, LWPOLY* poly, double tol, int* nfaces) { @@ -6005,3 +6072,913 @@ lwt_AddPolygon(LWT_TOPOLOGY* topo, LWPOLY* poly, double tol, int* nfaces) return ids; } + +/* + *---- polygonizer + */ + +/* An array of pointers to EDGERING structures */ +typedef struct LWT_ISO_EDGE_TABLE_T { + LWT_ISO_EDGE *edges; + int size; +} LWT_ISO_EDGE_TABLE; + +static int +compare_iso_edges_by_id(const void *si1, const void *si2) +{ + int a = ((LWT_ISO_EDGE *)si1)->edge_id; + int b = ((LWT_ISO_EDGE *)si2)->edge_id; + if ( a < b ) + return -1; + else if ( a > b ) + return 1; + else + return 0; +} + +static LWT_ISO_EDGE * +_lwt_getIsoEdgeById(LWT_ISO_EDGE_TABLE *tab, LWT_ELEMID id) +{ + LWT_ISO_EDGE key; + key.edge_id = id; + + void *match = bsearch( &key, tab->edges, tab->size, + sizeof(LWT_ISO_EDGE), + compare_iso_edges_by_id); + return match; +} + +typedef struct LWT_EDGERING_ELEM_T { + /* externally owned */ + LWT_ISO_EDGE *edge; + /* 0 if false, 1 if true */ + int left; +} LWT_EDGERING_ELEM; + +/* A ring of edges */ +typedef struct LWT_EDGERING_T { + /* Signed edge identifiers + * positive ones are walked in their direction, negative ones + * in the opposite direction */ + LWT_EDGERING_ELEM **elems; + /* Number of edges in the ring */ + int size; + int capacity; + /* Bounding box of the ring */ + GBOX *env; + /* Bounding box of the ring in GEOS format (for STRTree) */ + GEOSGeometry *genv; +} LWT_EDGERING; + +#define LWT_EDGERING_INIT(a) { \ + (a)->size = 0; \ + (a)->capacity = 1; \ + (a)->elems = lwalloc(sizeof(LWT_EDGERING_ELEM *) * (a)->capacity); \ + (a)->env = NULL; \ + (a)->genv = NULL; \ +} + +#define LWT_EDGERING_PUSH(a, r) { \ + if ( (a)->size + 1 > (a)->capacity ) { \ + (a)->capacity *= 2; \ + (a)->elems = lwrealloc((a)->elems, sizeof(LWT_EDGERING_ELEM *) * (a)->capacity); \ + } \ + /* lwdebug(1, "adding elem %d (%p) of edgering %p", (a)->size, (r), (a)); */ \ + (a)->elems[(a)->size++] = (r); \ +} + +#define LWT_EDGERING_CLEAN(a) { \ + int i; for (i=0; i<(a)->size; ++i) { \ + if ( (a)->elems[i] ) { \ + /* lwdebug(1, "freeing elem %d (%p) of edgering %p", i, (a)->elems[i], (a)); */ \ + lwfree((a)->elems[i]); \ + } \ + } \ + if ( (a)->elems ) { lwfree((a)->elems); (a)->elems = NULL; } \ + (a)->size = 0; \ + (a)->capacity = 0; \ + if ( (a)->env ) { lwfree((a)->env); (a)->env = NULL; } \ + if ( (a)->genv ) { GEOSGeom_destroy((a)->genv); (a)->genv = NULL; } \ +} + +/* An array of pointers to EDGERING structures */ +typedef struct LWT_EDGERING_ARRAY_T { + LWT_EDGERING **rings; + int size; + int capacity; + GEOSSTRtree* tree; +} LWT_EDGERING_ARRAY; + +#define LWT_EDGERING_ARRAY_INIT(a) { \ + (a)->size = 0; \ + (a)->capacity = 1; \ + (a)->rings = lwalloc(sizeof(LWT_EDGERING *) * (a)->capacity); \ + (a)->tree = NULL; \ +} + +/* WARNING: use of 'j' is intentional, not to clash with + * 'i' used in LWT_EDGERING_CLEAN */ +#define LWT_EDGERING_ARRAY_CLEAN(a) { \ + int j; for (j=0; j<(a)->size; ++j) { \ + LWT_EDGERING_CLEAN((a)->rings[j]); \ + } \ + if ( (a)->capacity ) lwfree((a)->rings); \ + if ( (a)->tree ) { \ + GEOSSTRtree_destroy( (a)->tree ); \ + (a)->tree = NULL; \ + } \ +} + +#define LWT_EDGERING_ARRAY_PUSH(a, r) { \ + if ( (a)->size + 1 > (a)->capacity ) { \ + (a)->capacity *= 2; \ + (a)->rings = lwrealloc((a)->rings, sizeof(LWT_EDGERING *) * (a)->capacity); \ + } \ + (a)->rings[(a)->size++] = (r); \ +} + +typedef struct LWT_EDGERING_POINT_ITERATOR_T { + LWT_EDGERING *ring; + LWT_EDGERING_ELEM *curelem; + int curelemidx; + int curidx; +} LWT_EDGERING_POINT_ITERATOR; + +static int +_lwt_EdgeRingIterator_next(LWT_EDGERING_POINT_ITERATOR *it, POINT2D *pt) +{ + LWT_EDGERING_ELEM *el = it->curelem; + POINTARRAY *pa; + + if ( ! el ) return 0; /* finished */ + + pa = el->edge->geom->points; + + int tonext = 0; + LWDEBUGF(3, "iterator fetching idx %d from pa of %d points", it->curidx, pa->npoints); + getPoint2d_p(pa, it->curidx, pt); + if ( el->left ) { + it->curidx++; + if ( it->curidx >= pa->npoints ) tonext = 1; + } else { + it->curidx--; + if ( it->curidx < 0 ) tonext = 1; + } + + if ( tonext ) + { + LWDEBUG(3, "iterator moving to next element"); + it->curelemidx++; + if ( it->curelemidx < it->ring->size ) + { + el = it->curelem = it->ring->elems[it->curelemidx]; + it->curidx = el->left ? 0 : el->edge->geom->points->npoints - 1; + } + else + { + it->curelem = NULL; + } + } + + return 1; +} + +/* Release return with lwfree */ +static LWT_EDGERING_POINT_ITERATOR * +_lwt_EdgeRingIterator_begin(LWT_EDGERING *er) +{ + LWT_EDGERING_POINT_ITERATOR *ret = lwalloc(sizeof(LWT_EDGERING_POINT_ITERATOR)); + ret->ring = er; + if ( er->size ) ret->curelem = er->elems[0]; + else ret->curelem = NULL; + ret->curelemidx = 0; + ret->curidx = ret->curelem->left ? 0 : ret->curelem->edge->geom->points->npoints - 1; + return ret; +} + +/* Identifier for a placeholder face that will be + * used to mark hole rings */ +#define LWT_HOLES_FACE_PLACEHOLDER INT32_MIN + +static int +_lwt_FetchNextUnvisitedEdge(LWT_TOPOLOGY *topo, LWT_ISO_EDGE_TABLE *etab, int from) +{ + while ( + from < etab->size && + etab->edges[from].face_left != -1 && + etab->edges[from].face_right != -1 + ) from++; + return from < etab->size ? from : -1; +} + +static LWT_ISO_EDGE * +_lwt_FetchAllEdges(LWT_TOPOLOGY *topo, int *numedges) +{ + LWT_ISO_EDGE *edge; + int fields = LWT_COL_EDGE_ALL; + int nelems = 1; + + edge = lwt_be_getEdgeWithinBox2D( topo, NULL, &nelems, fields, 0); + *numedges = nelems; + if ( nelems == -1 ) { + lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); + return NULL; + } + return edge; +} + +/* Update the side face of given ring edges + * + * Edge identifiers are signed, those with negative identifier + * need to be updated their right_face, those with positive + * identifier need to be updated their left_face. + * + * @param face identifier of the face bound by the ring + * @return 0 on success, -1 on error + */ +static int +_lwt_UpdateEdgeRingSideFace(LWT_TOPOLOGY *topo, LWT_EDGERING *ring, + LWT_ELEMID face) +{ + LWT_ISO_EDGE *forward_edges = NULL; + int forward_edges_count = 0; + LWT_ISO_EDGE *backward_edges = NULL; + int backward_edges_count = 0; + int i, ret; + + /* Make a list of forward_edges and backward_edges */ + + forward_edges = lwalloc(sizeof(LWT_ISO_EDGE) * ring->size); + forward_edges_count = 0; + backward_edges = lwalloc(sizeof(LWT_ISO_EDGE) * ring->size); + backward_edges_count = 0; + + for ( i=0; isize; ++i ) + { + LWT_EDGERING_ELEM *elem = ring->elems[i]; + LWT_ISO_EDGE *edge = elem->edge; + LWT_ELEMID id = edge->edge_id; + if ( elem->left ) + { + LWDEBUGF(3, "Forward edge %d is %d", forward_edges_count, id); + forward_edges[forward_edges_count].edge_id = id; + forward_edges[forward_edges_count++].face_left = face; + edge->face_left = face; + } + else + { + LWDEBUGF(3, "Backward edge %d is %d", forward_edges_count, id); + backward_edges[backward_edges_count].edge_id = id; + backward_edges[backward_edges_count++].face_right = face; + edge->face_right = face; + } + } + + /* Update forward edges */ + if ( forward_edges_count ) + { + ret = lwt_be_updateEdgesById(topo, forward_edges, + forward_edges_count, + LWT_COL_EDGE_FACE_LEFT); + if ( ret == -1 ) + { + lwfree( forward_edges ); + lwfree( backward_edges ); + lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); + return -1; + } + if ( ret != forward_edges_count ) + { + lwfree( forward_edges ); + lwfree( backward_edges ); + lwerror("Unexpected error: %d edges updated when expecting %d (forward)", + ret, forward_edges_count); + return -1; + } + } + + /* Update backward edges */ + if ( backward_edges_count ) + { + ret = lwt_be_updateEdgesById(topo, backward_edges, + backward_edges_count, + LWT_COL_EDGE_FACE_RIGHT); + if ( ret == -1 ) + { + lwfree( forward_edges ); + lwfree( backward_edges ); + lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); + return -1; + } + if ( ret != backward_edges_count ) + { + lwfree( forward_edges ); + lwfree( backward_edges ); + lwerror("Unexpected error: %d edges updated when expecting %d (backward)", + ret, backward_edges_count); + return -1; + } + } + + lwfree( forward_edges ); + lwfree( backward_edges ); + + return 0; +} + +/* + * @param side 1 for left side, -1 for right side + */ +static LWT_EDGERING * +_lwt_BuildEdgeRing(LWT_TOPOLOGY *topo, LWT_ISO_EDGE_TABLE *edges, + LWT_ISO_EDGE *edge, int side) +{ + LWT_EDGERING *ring; + LWT_EDGERING_ELEM *elem; + LWT_ISO_EDGE *cur; + int curside; + + ring = lwalloc(sizeof(LWT_EDGERING)); + LWT_EDGERING_INIT(ring); + + cur = edge; + curside = side; + + LWDEBUGF(2, "Building rings for edge %d (side %d)", cur->edge_id, side); + + do { + LWT_ELEMID next; + + elem = lwalloc(sizeof(LWT_EDGERING_ELEM)); + elem->edge = cur; + elem->left = ( curside == 1 ); + + /* Mark edge as "visited" */ + if ( elem->left ) cur->face_left = LWT_HOLES_FACE_PLACEHOLDER; + else cur->face_right = LWT_HOLES_FACE_PLACEHOLDER; + + LWT_EDGERING_PUSH(ring, elem); + next = elem->left ? cur->next_left : cur->next_right; + + LWDEBUGF(3, " next edge is %d", next); + + if ( next > 0 ) curside = 1; + else { curside = -1; next = -next; } + cur = _lwt_getIsoEdgeById(edges, next); + if ( ! cur ) + { + lwerror("Could not find edge with id %d", next); + break; + } + } while (cur != edge || curside != side); + + LWDEBUGF(1, "Ring for edge %d has %d elems", edge->edge_id*side, ring->size); + + return ring; +} + +static double +_lwt_EdgeRingSignedArea(LWT_EDGERING_POINT_ITERATOR *it) +{ + POINT2D P1; + POINT2D P2; + POINT2D P3; + double sum = 0.0; + double x0, x, y1, y2; + + if ( ! _lwt_EdgeRingIterator_next(it, &P1) ) return 0.0; + if ( ! _lwt_EdgeRingIterator_next(it, &P2) ) return 0.0; + + LWDEBUG(2, "_lwt_EdgeRingSignedArea"); + + x0 = P1.x; + while ( _lwt_EdgeRingIterator_next(it, &P3) ) + { + x = P2.x - x0; + y1 = P3.y; + y2 = P1.y; + sum += x * (y2-y1); + + /* Move forwards! */ + P1 = P2; + P2 = P3; + } + + return sum / 2.0; +} + + +/* Return 1 for true, 0 for false */ +static int +_lwt_EdgeRingIsCCW(LWT_EDGERING *ring) +{ + double sa; + + LWDEBUGF(2, "_lwt_EdgeRingIsCCW, ring has %d elems", ring->size); + LWT_EDGERING_POINT_ITERATOR *it = _lwt_EdgeRingIterator_begin(ring); + sa = _lwt_EdgeRingSignedArea(it); + LWDEBUGF(2, "_lwt_EdgeRingIsCCW, signed area is %g", sa); + lwfree(it); + if ( sa >= 0 ) return 0; + else return 1; +} + +static int +_lwt_EdgeRingCrossingCount(const POINT2D *p, LWT_EDGERING_POINT_ITERATOR *it) +{ + int cn = 0; /* the crossing number counter */ + POINT2D v1, v2; +#ifndef RELAX + POINT2D v0; +#endif + + if ( ! _lwt_EdgeRingIterator_next(it, &v1) ) return cn; + v0 = v1; + while ( _lwt_EdgeRingIterator_next(it, &v2) ) + { + double vt; + + /* edge from vertex i to vertex i+1 */ + if + ( + /* an upward crossing */ + ((v1.y <= p->y) && (v2.y > p->y)) + /* a downward crossing */ + || ((v1.y > p->y) && (v2.y <= p->y)) + ) + { + + vt = (double)(p->y - v1.y) / (v2.y - v1.y); + + /* P->x x < v1.x + vt * (v2.x - v1.x)) + { + /* a valid crossing of y=p->y right of p->x */ + ++cn; + } + } + v1 = v2; + } + + LWDEBUGF(3, "_lwt_EdgeRingCrossingCount returning %d", cn); + +#ifndef RELAX + if ( memcmp(&v1, &v0, sizeof(POINT2D)) ) + { + lwerror("_lwt_EdgeRingCrossingCount: V[n] != V[0] (%g %g != %g %g)", + v1.x, v1.y, v0.x, v0.y); + return -1; + } +#endif + + return cn; +} + +/* Return 1 for true, 0 for false */ +static int +_lwt_EdgeRingContainsPoint(LWT_EDGERING *ring, POINT2D *p) +{ + int cn = 0; + + LWT_EDGERING_POINT_ITERATOR *it = _lwt_EdgeRingIterator_begin(ring); + cn = _lwt_EdgeRingCrossingCount(p, it); + lwfree(it); + return (cn&1); /* 0 if even (out), and 1 if odd (in) */ +} + +static GBOX * +_lwt_EdgeRingGetBbox(LWT_EDGERING *ring) +{ + int i; + + if ( ! ring->env ) + { + LWDEBUGF(2, "Computing GBOX for ring %p", ring); + for (i=0; isize; ++i) + { + LWT_EDGERING_ELEM *elem = ring->elems[i]; + LWLINE *g = elem->edge->geom; + const GBOX *newbox = lwgeom_get_bbox(lwline_as_lwgeom(g)); + if ( ! i ) ring->env = gbox_clone( newbox ); + else gbox_merge( newbox, ring->env ); + } + } + + return ring->env; +} + +static LWT_ELEMID +_lwt_EdgeRingGetFace(LWT_EDGERING *ring) +{ + LWT_EDGERING_ELEM *el = ring->elems[0]; + return el->left ? el->edge->face_left : el->edge->face_right; +} + + +/* + * Register a face on an edge side + * + * Create and register face to shell (CCW) walks, + * register arbitrary negative face_id to CW rings. + * + * Push CCW rings to shells, CW rings to holes. + * + * The ownership of the "geom" and "ids" members of the + * LWT_EDGERING pushed to the given LWT_EDGERING_ARRAYS + * are transferred to caller. + * + * @param side 1 for left side, -1 for right side + * + * @param holes an array where holes will be pushed + * + * @param shells an array where shells will be pushed + * + * @param registered id of registered face. It will be a negative number + * for holes or isolated edge strips (still registered in the face + * table, but only temporary). + * + * @return 0 on success, -1 on error. + * + */ +static int +_lwt_RegisterFaceOnEdgeSide(LWT_TOPOLOGY *topo, LWT_ISO_EDGE *edge, + int side, LWT_ISO_EDGE_TABLE *edges, + LWT_EDGERING_ARRAY *holes, + LWT_EDGERING_ARRAY *shells, + LWT_ELEMID *registered) +{ + const LWT_BE_IFACE *iface = topo->be_iface; + int sedge = edge->edge_id * side; + /* this is arbitrary, could be taken as parameter */ + static const int placeholder_faceid = LWT_HOLES_FACE_PLACEHOLDER; + LWT_EDGERING *ring; + + /* Get edge ring */ + ring = _lwt_BuildEdgeRing(topo, edges, edge, side); + + LWDEBUG(2, "Ring built, calling EdgeRingIsCCW"); + + /* Compute winding (CW or CCW?) */ + int isccw = _lwt_EdgeRingIsCCW(ring); + + if ( isccw ) + { + /* Create new face */ + LWT_ISO_FACE newface; + + LWDEBUGF(1, "Ring of edge %d is a shell (shell %d)", sedge, shells->size); + + newface.mbr = _lwt_EdgeRingGetBbox(ring); + + newface.face_id = -1; + /* Insert the new face */ + int ret = lwt_be_insertFaces( topo, &newface, 1 ); + newface.mbr = NULL; + if ( ret == -1 ) + { + lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); + return -1; + } + if ( ret != 1 ) + { + lwerror("Unexpected error: %d faces inserted when expecting 1", ret); + return -1; + } + /* return new face_id */ + *registered = newface.face_id; + LWT_EDGERING_ARRAY_PUSH(shells, ring); + + /* update ring edges set new face_id on resp. side to *registered */ + ret = _lwt_UpdateEdgeRingSideFace(topo, ring, *registered); + if ( ret ) + { + lwerror("Errors updating edgering side face: %s", + lwt_be_lastErrorMessage(iface)); + return -1; + } + + } + else /* cw, so is an hole */ + { + LWDEBUGF(1, "Ring of edge %d is a hole (hole %d)", sedge, holes->size); + *registered = placeholder_faceid; + LWT_EDGERING_ARRAY_PUSH(holes, ring); + } + + return 0; +} + +static void +_lwt_AccumulateCanditates(void* item, void* userdata) +{ + LWT_EDGERING_ARRAY *candidates = userdata; + LWT_EDGERING *sring = item; + LWT_EDGERING_ARRAY_PUSH(candidates, sring); +} + +static LWT_ELEMID +_lwt_FindFaceContainingRing(LWT_TOPOLOGY* topo, LWT_EDGERING *ring, + LWT_EDGERING_ARRAY *shells) +{ + LWT_ELEMID foundInFace = -1; + int i; + const GBOX *minenv = NULL; + POINT2D pt; + const GBOX *testbox; + GEOSGeometry *ghole; + + getPoint2d_p( ring->elems[0]->edge->geom->points, 0, &pt ); + + testbox = _lwt_EdgeRingGetBbox(ring); + + /* Create a GEOS Point from a vertex of the hole ring */ + { + LWPOINT *point = lwpoint_make2d(topo->srid, pt.x, pt.y); + ghole = LWGEOM2GEOS( lwpoint_as_lwgeom(point), 1 ); + lwpoint_free(point); + if ( ! ghole ) { + lwerror("Could not convert edge geometry to GEOS: %s", lwgeom_geos_errmsg); + return -1; + } + } + + /* Build STRtree of shell envelopes */ + if ( ! shells->tree ) + { + static const int STRTREE_NODE_CAPACITY = 10; + LWDEBUG(1, "Building STRtree"); + shells->tree = GEOSSTRtree_create(STRTREE_NODE_CAPACITY); + if (shells->tree == NULL) + { + lwerror("Could not create GEOS STRTree: %s", lwgeom_geos_errmsg); + return -1; + } + for (i=0; isize; ++i) + { + LWT_EDGERING *sring = shells->rings[i]; + const GBOX* shellbox = _lwt_EdgeRingGetBbox(sring); + LWDEBUGF(2, "GBOX of shell %p for edge %d is %g %g,%g %g", + sring, sring->elems[0]->edge->edge_id, shellbox->xmin, + shellbox->ymin, shellbox->xmax, shellbox->ymax); + POINTARRAY *pa = ptarray_construct(0, 0, 2); + POINT4D pt; + LWLINE *diag; + pt.x = shellbox->xmin; + pt.y = shellbox->ymin; + ptarray_set_point4d(pa, 0, &pt); + pt.x = shellbox->xmax; + pt.y = shellbox->ymax; + ptarray_set_point4d(pa, 1, &pt); + diag = lwline_construct(topo->srid, NULL, pa); + /* Record just envelope in ggeom */ + /* making valid, probably not needed */ + sring->genv = LWGEOM2GEOS( lwline_as_lwgeom(diag), 1 ); + lwline_free(diag); + GEOSSTRtree_insert(shells->tree, sring->genv, sring); + } + LWDEBUG(1, "STRtree build completed"); + } + + LWT_EDGERING_ARRAY candidates; + LWT_EDGERING_ARRAY_INIT(&candidates); + GEOSSTRtree_query(shells->tree, ghole, &_lwt_AccumulateCanditates, &candidates); + LWDEBUGF(1, "Found %d candidate shells containing first point of ring's originating edge %d", + candidates.size, ring->elems[0]->edge->edge_id * ( ring->elems[0]->left ? 1 : -1 ) ); + + /* TODO: sort candidates by bounding box size */ + + for (i=0; ielems[0]->edge->edge_id == ring->elems[0]->edge->edge_id ) + { + LWDEBUGF(1, "Shell %d is on other side of ring", + _lwt_EdgeRingGetFace(sring)); + continue; + } + + /* The hole envelope cannot equal the shell envelope */ + if ( gbox_same(shellbox, testbox) ) + { + LWDEBUGF(1, "Bbox of shell %d equals that of hole ring", + _lwt_EdgeRingGetFace(sring)); + continue; + } + + /* Skip if ring box is not in shell box */ + if ( ! gbox_contains_2d(shellbox, testbox) ) + { + LWDEBUGF(1, "Bbox of shell %d does not contain bbox of ring point", + _lwt_EdgeRingGetFace(sring)); + continue; + } + + /* Skip test if a containing shell was already found + * and this shell's bbox is not contained in the other */ + if ( minenv && ! gbox_contains_2d(minenv, shellbox) ) + { + LWDEBUGF(2, "Bbox of shell %d (face %d) not contained by bbox " + "of last shell found to contain the point", + i, _lwt_EdgeRingGetFace(sring)); + continue; + } + + contains = _lwt_EdgeRingContainsPoint(sring, &pt); + if ( contains ) + { + /* Continue until all shells are tested, as we want to + * use the one with the smallest bounding box */ + /* IDEA: sort shells by bbox size, stopping on first match */ + LWDEBUGF(1, "Shell %d contains hole of edge %d", + _lwt_EdgeRingGetFace(sring), + ring->elems[0]->edge->edge_id); + minenv = shellbox; + foundInFace = _lwt_EdgeRingGetFace(sring); + } + } + if ( foundInFace == -1 ) foundInFace = 0; + + candidates.size = 0; /* Avoid destroying the actual shell rings */ + LWT_EDGERING_ARRAY_CLEAN(&candidates); + + GEOSGeom_destroy(ghole); + + return foundInFace; +} + +/* + * @return -1 on error (and report error), + * 1 if faces beside the universal one exist + * 0 otherwise + */ +static int +_lwt_CheckFacesExist(LWT_TOPOLOGY *topo) +{ + LWT_ISO_FACE *faces; + int fields = LWT_COL_FACE_FACE_ID; + int nelems = 1; + GBOX qbox; + + qbox.xmin = qbox.ymin = -DBL_MAX; + qbox.xmax = qbox.ymax = DBL_MAX; + faces = lwt_be_getFaceWithinBox2D( topo, &qbox, &nelems, fields, 1); + if ( nelems == -1 ) { + lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); + return -1; + } + if ( faces ) _lwt_release_faces(faces, nelems); + return nelems; +} + +int +lwt_Polygonize(LWT_TOPOLOGY* topo) +{ + /* + Fetch all edges + Sort edges by edge_id + Mark all edges' left and right face as -1 + Iteratively: + Fetch next edge with left or right face == -1 + For each side with face == -1: + Find ring on its side + If ring is CCW: + create a new face, assign to the ring edges' appropriate side + If ring is CW (face needs to be same of external): + assign a negative face_id the ring edges' appropriate side + Now for each edge with a negative face_id on the side: + Find containing face (mbr cache and all) + Update with id of containing face + */ + + const LWT_BE_IFACE *iface = topo->be_iface; + LWT_ISO_EDGE *edge; + int numfaces = -1; + LWT_ISO_EDGE_TABLE edgetable; + LWT_EDGERING_ARRAY holes, shells; + int i; + int err = 0; + + LWT_EDGERING_ARRAY_INIT(&holes); + LWT_EDGERING_ARRAY_INIT(&shells); + + initGEOS(lwnotice, lwgeom_geos_error); + + /* + Check if Topology already contains some Face + (ignoring the Universal Face) + */ + numfaces = _lwt_CheckFacesExist(topo); + if ( numfaces != 0 ) { + if ( numfaces > 0 ) { + /* Faces exist */ + lwerror("Polygonize: face table is not empty."); + } + /* Backend error, message should have been printed already */ + return -1; + } + + + edgetable.edges = _lwt_FetchAllEdges(topo, &(edgetable.size)); + if ( ! edgetable.edges ) { + if (edgetable.size == 0) { + /* not an error: no Edges */ + return 0; + } + /* error should have been printed already */ + return -1; + } + + /* Sort edges by ID (to allow btree searches) */ + qsort(edgetable.edges, edgetable.size, sizeof(LWT_ISO_EDGE), compare_iso_edges_by_id); + + /* Mark all edges as unvisited */ + for (i=0; iedge_id, edge->face_left, edge->face_right); + if ( edge->face_left == -1 ) + { + err = _lwt_RegisterFaceOnEdgeSide(topo, edge, 1, &edgetable, + &holes, &shells, &newface); + if ( err ) break; + LWDEBUGF(1, "New face on the left of edge %d is %d", + edge->edge_id, newface); + edge->face_left = newface; + } + if ( edge->face_right == -1 ) + { + err = _lwt_RegisterFaceOnEdgeSide(topo, edge, -1, &edgetable, + &holes, &shells, &newface); + if ( err ) break; + LWDEBUGF(1, "New face on the right of edge %d is %d", + edge->edge_id, newface); + edge->face_right = newface; + } + } + + if ( err ) + { + _lwt_release_edges(edgetable.edges, edgetable.size); + LWT_EDGERING_ARRAY_CLEAN( &holes ); + LWT_EDGERING_ARRAY_CLEAN( &shells ); + lwerror("Errors fetching or registering face-missing edges: %s", + lwt_be_lastErrorMessage(iface)); + return -1; + } + + LWDEBUGF(1, "Found %d holes and %d shells", holes.size, shells.size); + + /* TODO: sort holes by pt.x, sort shells by bbox.xmin */ + + /* Assign shells to holes */ + for (i=0; i