[postgis-tickets] r15857 - Add lwt_Polygonize and lwt_AddLineNoFace functions
Sandro Santilli
strk at kbt.io
Sat Sep 30 08:26:28 PDT 2017
Author: strk
Date: 2017-09-30 08:26:28 -0700 (Sat, 30 Sep 2017)
New Revision: 15857
Modified:
trunk/liblwgeom/liblwgeom_topo.h
trunk/liblwgeom/lwgeom_topo.c
Log:
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.
Modified: trunk/liblwgeom/liblwgeom_topo.h
===================================================================
--- trunk/liblwgeom/liblwgeom_topo.h 2017-09-30 06:38:52 UTC (rev 15856)
+++ trunk/liblwgeom/liblwgeom_topo.h 2017-09-30 15:26:28 UTC (rev 15857)
@@ -18,7 +18,7 @@
*
**********************************************************************
*
- * Copyright (C) 2015 Sandro Santilli <strk at kbt.io>
+ * Copyright (C) 2015-2017 Sandro Santilli <strk at kbt.io>
*
**********************************************************************/
@@ -1057,6 +1057,60 @@
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 <nedges> 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
*
* The boundary of the given polygon will snap to existing nodes or
Modified: trunk/liblwgeom/lwgeom_topo.c
===================================================================
--- trunk/liblwgeom/lwgeom_topo.c 2017-09-30 06:38:52 UTC (rev 15856)
+++ trunk/liblwgeom/lwgeom_topo.c 2017-09-30 15:26:28 UTC (rev 15857)
@@ -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 @@
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 @@
}
}
- if ( face == -1 || ! skipISOChecks )
+ if ( checkFace && ( face == -1 || ! skipISOChecks ) )
{
foundInFace = lwt_be_getFaceContainingPoint(topo, pt); /*x*/
if ( foundInFace == -2 ) {
@@ -577,6 +585,13 @@
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,62 +1769,13 @@
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; i<num_signed_edge_ids; ++i) {
- if ( signed_edge_ids[i] == -sedge ) {
- /* No split here */
- LWDEBUG(1, "not a ring");
- lwfree( signed_edge_ids );
- return 0;
- }
- }
-
- LWDEBUGF(1, "Edge %" LWTFMT_ELEMID " split face %" LWTFMT_ELEMID " (mbr_only:%d)",
- sedge, face, mbr_only);
-
/* Construct a polygon using edges of the ring */
numedges = 0;
edge_ids = lwalloc(sizeof(LWT_ELEMID)*num_signed_edge_ids);
@@ -1831,17 +1797,15 @@
lwfree( edge_ids );
if ( i == -1 )
{
- lwfree( signed_edge_ids );
- /* ring_edges should be NULL */
lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_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 @@
for ( i=0; i<num_signed_edge_ids; ++i )
{
LWT_ELEMID eid = signed_edge_ids[i];
- LWDEBUGF(1, "Edge %d in ring of edge %" LWTFMT_ELEMID " is edge %" LWTFMT_ELEMID,
- i, sedge, eid);
+ LWDEBUGF(2, "Edge %d in ring is edge %" LWTFMT_ELEMID, i, eid);
LWT_ISO_EDGE *edge = NULL;
POINTARRAY *epa;
for ( j=0; j<numedges; ++j )
@@ -1865,10 +1828,9 @@
}
if ( edge == NULL )
{
- lwfree( signed_edge_ids );
_lwt_release_edges(ring_edges, numedges);
lwerror("missing edge that was found in ring edges loop");
- return -2;
+ return NULL;
}
if ( pa == NULL )
@@ -1892,13 +1854,82 @@
}
}
}
+ _lwt_release_edges(ring_edges, numedges);
POINTARRAY **points = lwalloc(sizeof(POINTARRAY*));
points[0] = pa;
+
/* NOTE: the ring may very well have collapsed components,
* which would make it topologically invalid
*/
LWPOLY* shell = lwpoly_construct(0, 0, 1, points);
+ return shell;
+}
+/*
+ * 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 )
+{
+ 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; i<num_signed_edge_ids; ++i) {
+ if ( signed_edge_ids[i] == -sedge ) {
+ /* No split here */
+ LWDEBUG(1, "not a ring");
+ lwfree( signed_edge_ids );
+ return 0;
+ }
+ }
+
+ LWDEBUGF(1, "Edge %" LWTFMT_ELEMID " split face %" LWTFMT_ELEMID " (mbr_only:%d)",
+ sedge, face, mbr_only);
+
+ /* Construct a polygon using edges of the ring */
+ LWPOLY *shell = _lwt_MakeRingShell(topo, signed_edge_ids,
+ num_signed_edge_ids);
+ if ( ! shell ) {
+ lwfree( signed_edge_ids );
+ /* ring_edges should be NULL */
+ lwerror("Could not create ring shell: %s", lwt_be_lastErrorMessage(topo->be_iface));
+ return -2;
+ }
+ const POINTARRAY *pa = shell->rings[0];
+
int isccw = ptarray_isccw(pa);
LWDEBUGF(1, "Ring of edge %" LWTFMT_ELEMID " is %sclockwise",
sedge, isccw ? "counter" : "");
@@ -1911,7 +1942,6 @@
{
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 @@
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 @@
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 @@
{
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 @@
{
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 @@
{
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 @@
{
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 @@
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 @@
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 @@
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 @@
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 @@
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 @@
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_release_edges(ring_edges, numedges);
_lwt_release_edges(edges, numfaceedges);
/* Update isolated nodes which are now in new face */
@@ -2335,6 +2351,13 @@
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 @@
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 @@
}
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 @@
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,8 +2743,10 @@
}
}
- /* Check face splitting */
+ /* Check face splitting, if required */
+ if ( modFace > -1 ) {
+
if ( ! isclosed && ( epan.was_isolated || span.was_isolated ) )
{
LWDEBUG(1, "New edge is dangling, so it cannot split any face");
@@ -2788,6 +2813,8 @@
}
}
+ } // end of face split checking
+
return newedge.edge_id;
}
@@ -3601,75 +3628,82 @@
-- 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 )
+ 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 left of edge %" LWTFMT_ELEMID,
- oldedge->face_left, edge_id);
- return -1;
+ 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 0
+ 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(nface1, WKT_EXTENDED, 2, &sz);
- LWDEBUGF(1, "new geometry of face left (%d): %s", (int)oldedge->face_left, 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(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 )
+ LWDEBUGF(1, "%d faces to update", facestoupdate);
+ if ( facestoupdate )
{
- lwerror("lwt_ChangeEdgeGeom could not construct face %"
- LWTFMT_ELEMID ", on the right of edge %" LWTFMT_ELEMID,
- oldedge->face_right, edge_id);
- 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 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 */
- }
- LWDEBUGF(1, "%d faces to update", facestoupdate);
- if ( facestoupdate )
- {
- 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 @@
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 @@
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 @@
{
/* 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 @@
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 @@
/*
* 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 @@
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 @@
"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 @@
/* 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 @@
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 @@
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_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 @@
}
#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 )
{
@@ -5872,6 +5927,18 @@
}
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)
{
int i;
@@ -6005,3 +6072,913 @@
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; i<ring->size; ++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 <intersect */
+ if (p->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; i<ring->size; ++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; i<shells->size; ++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; i<candidates.size; ++i)
+ {
+ LWT_EDGERING *sring = candidates.rings[i];
+ const GBOX* shellbox = _lwt_EdgeRingGetBbox(sring);
+ int contains = 0;
+
+ if ( sring->elems[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; i<edgetable.size; ++i)
+ edgetable.edges[i].face_left = edgetable.edges[i].face_right = -1;
+
+ i = 0;
+ while (1)
+ {
+ i = _lwt_FetchNextUnvisitedEdge(topo, &edgetable, i);
+ if ( i < 0 ) break; /* end of unvisited */
+ edge = &(edgetable.edges[i]);
+
+ LWT_ELEMID newface = -1;
+
+ LWDEBUGF(1, "Next face-missing edge has id:%d, face_left:%d, face_right:%d",
+ edge->edge_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<holes.size; ++i)
+ {
+ LWT_ELEMID containing_face;
+ LWT_EDGERING *ring = holes.rings[i];
+
+ containing_face = _lwt_FindFaceContainingRing(topo, ring, &shells);
+ LWDEBUGF(1, "Ring %d contained by face %" LWTFMT_ELEMID, i, containing_face);
+ if ( containing_face == -1 )
+ {
+ _lwt_release_edges(edgetable.edges, edgetable.size);
+ LWT_EDGERING_ARRAY_CLEAN( &holes );
+ LWT_EDGERING_ARRAY_CLEAN( &shells );
+ lwerror("Errors finding face containing ring: %s",
+ lwt_be_lastErrorMessage(iface));
+ return -1;
+ }
+ int ret = _lwt_UpdateEdgeRingSideFace(topo, holes.rings[i], containing_face);
+ if ( ret )
+ {
+ _lwt_release_edges(edgetable.edges, edgetable.size);
+ LWT_EDGERING_ARRAY_CLEAN( &holes );
+ LWT_EDGERING_ARRAY_CLEAN( &shells );
+ lwerror("Errors updating edgering side face: %s",
+ lwt_be_lastErrorMessage(iface));
+ return -1;
+ }
+ }
+
+ LWDEBUG(1, "All holes assigned, cleaning up");
+
+ _lwt_release_edges(edgetable.edges, edgetable.size);
+
+ /* delete all shell and hole EDGERINGS */
+ LWT_EDGERING_ARRAY_CLEAN( &holes );
+ LWT_EDGERING_ARRAY_CLEAN( &shells );
+
+ return 0;
+}
More information about the postgis-tickets
mailing list