/********************************************************************** * * PostGIS - Spatial Types for PostgreSQL * http://postgis.net * * PostGIS is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 2 of the License, or * (at your option) any later version. * * PostGIS is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with PostGIS. If not, see . * ********************************************************************** * * Copyright (C) 2015-2022 Sandro Santilli * **********************************************************************/ #include "../postgis_config.h" /*#define POSTGIS_DEBUG_LEVEL 1*/ #include "lwgeom_log.h" #include "liblwgeom_internal.h" #include "liblwgeom_topo_internal.h" #include "lwgeom_geos.h" #include #include /* for PRId64 */ #include #ifdef WIN32 # define LWTFMT_ELEMID "lld" #else # define LWTFMT_ELEMID PRId64 #endif /********************************************************************* * * Backend iface * ********************************************************************/ LWT_BE_IFACE* lwt_CreateBackendIface(const LWT_BE_DATA *data) { LWT_BE_IFACE *iface = lwalloc(sizeof(LWT_BE_IFACE)); iface->data = data; iface->cb = NULL; return iface; } void lwt_BackendIfaceRegisterCallbacks(LWT_BE_IFACE *iface, const LWT_BE_CALLBACKS* cb) { iface->cb = cb; } void lwt_FreeBackendIface(LWT_BE_IFACE* iface) { lwfree(iface); } /********************************************************************* * * Backend wrappers * ********************************************************************/ #define CHECKCB(be, method) do { \ if ( ! (be)->cb || ! (be)->cb->method ) \ lwerror("Callback " # method " not registered by backend"); \ } while (0) #define CB0(be, method) \ CHECKCB(be, method);\ return (be)->cb->method((be)->data) #define CB1(be, method, a1) \ CHECKCB(be, method);\ return (be)->cb->method((be)->data, a1) #define CBT0(to, method) \ CHECKCB((to)->be_iface, method);\ return (to)->be_iface->cb->method((to)->be_topo) #define CBT1(to, method, a1) \ CHECKCB((to)->be_iface, method);\ return (to)->be_iface->cb->method((to)->be_topo, a1) #define CBT2(to, method, a1, a2) \ CHECKCB((to)->be_iface, method);\ return (to)->be_iface->cb->method((to)->be_topo, a1, a2) #define CBT3(to, method, a1, a2, a3) \ CHECKCB((to)->be_iface, method);\ return (to)->be_iface->cb->method((to)->be_topo, a1, a2, a3) #define CBT4(to, method, a1, a2, a3, a4) \ CHECKCB((to)->be_iface, method);\ return (to)->be_iface->cb->method((to)->be_topo, a1, a2, a3, a4) #define CBT5(to, method, a1, a2, a3, a4, a5) \ CHECKCB((to)->be_iface, method);\ return (to)->be_iface->cb->method((to)->be_topo, a1, a2, a3, a4, a5) #define CBT6(to, method, a1, a2, a3, a4, a5, a6) \ CHECKCB((to)->be_iface, method);\ return (to)->be_iface->cb->method((to)->be_topo, a1, a2, a3, a4, a5, a6) const char * lwt_be_lastErrorMessage(const LWT_BE_IFACE* be) { CB0(be, lastErrorMessage); } LWT_BE_TOPOLOGY * lwt_be_loadTopologyByName(LWT_BE_IFACE *be, const char *name) { CB1(be, loadTopologyByName, name); } static int lwt_be_topoGetSRID(LWT_TOPOLOGY *topo) { CBT0(topo, topoGetSRID); } static double lwt_be_topoGetPrecision(LWT_TOPOLOGY *topo) { CBT0(topo, topoGetPrecision); } static int lwt_be_topoHasZ(LWT_TOPOLOGY *topo) { CBT0(topo, topoHasZ); } int lwt_be_freeTopology(LWT_TOPOLOGY *topo) { CBT0(topo, freeTopology); } LWT_ISO_NODE * lwt_be_getNodeById(LWT_TOPOLOGY *topo, const LWT_ELEMID *ids, uint64_t *numelems, int fields) { CBT3(topo, getNodeById, ids, numelems, fields); } LWT_ISO_NODE * lwt_be_getNodeWithinDistance2D(LWT_TOPOLOGY *topo, LWPOINT *pt, double dist, uint64_t *numelems, int fields, int64_t limit) { CBT5(topo, getNodeWithinDistance2D, pt, dist, numelems, fields, limit); } static LWT_ISO_NODE * lwt_be_getNodeWithinBox2D(const LWT_TOPOLOGY *topo, const GBOX *box, uint64_t *numelems, int fields, uint64_t limit) { CBT4(topo, getNodeWithinBox2D, box, numelems, fields, limit); } static LWT_ISO_EDGE * lwt_be_getEdgeWithinBox2D(const LWT_TOPOLOGY *topo, const GBOX *box, uint64_t *numelems, int fields, uint64_t limit) { CBT4(topo, getEdgeWithinBox2D, box, numelems, fields, limit); } static LWT_ISO_FACE * lwt_be_getFaceWithinBox2D(const LWT_TOPOLOGY *topo, const GBOX *box, uint64_t *numelems, int fields, uint64_t limit) { CBT4(topo, getFaceWithinBox2D, box, numelems, fields, limit); } int lwt_be_insertNodes(LWT_TOPOLOGY *topo, LWT_ISO_NODE *node, uint64_t numelems) { CBT2(topo, insertNodes, node, numelems); } static int lwt_be_insertFaces(LWT_TOPOLOGY *topo, LWT_ISO_FACE *face, uint64_t numelems) { CBT2(topo, insertFaces, face, numelems); } static int lwt_be_deleteFacesById(const LWT_TOPOLOGY *topo, const LWT_ELEMID *ids, uint64_t numelems) { CBT2(topo, deleteFacesById, ids, numelems); } static int lwt_be_deleteNodesById(const LWT_TOPOLOGY *topo, const LWT_ELEMID *ids, uint64_t numelems) { CBT2(topo, deleteNodesById, ids, numelems); } LWT_ELEMID lwt_be_getNextEdgeId(LWT_TOPOLOGY* topo) { CBT0(topo, getNextEdgeId); } LWT_ISO_EDGE * lwt_be_getEdgeById(LWT_TOPOLOGY *topo, const LWT_ELEMID *ids, uint64_t *numelems, int fields) { CBT3(topo, getEdgeById, ids, numelems, fields); } static LWT_ISO_FACE * lwt_be_getFaceById(LWT_TOPOLOGY *topo, const LWT_ELEMID *ids, uint64_t *numelems, int fields) { CBT3(topo, getFaceById, ids, numelems, fields); } static LWT_ISO_EDGE * lwt_be_getEdgeByNode(LWT_TOPOLOGY *topo, const LWT_ELEMID *ids, uint64_t *numelems, int fields) { CBT3(topo, getEdgeByNode, ids, numelems, fields); } static LWT_ISO_EDGE * lwt_be_getEdgeByFace(LWT_TOPOLOGY *topo, const LWT_ELEMID *ids, uint64_t *numelems, int fields, const GBOX *box) { CBT4(topo, getEdgeByFace, ids, numelems, fields, box); } static LWT_ISO_NODE * lwt_be_getNodeByFace(LWT_TOPOLOGY *topo, const LWT_ELEMID *ids, uint64_t *numelems, int fields, const GBOX *box) { CBT4(topo, getNodeByFace, ids, numelems, fields, box); } LWT_ISO_EDGE * lwt_be_getEdgeWithinDistance2D(LWT_TOPOLOGY *topo, const LWPOINT *pt, double dist, uint64_t *numelems, int fields, int64_t limit) { CBT5(topo, getEdgeWithinDistance2D, pt, dist, numelems, fields, limit); } int lwt_be_insertEdges(LWT_TOPOLOGY *topo, LWT_ISO_EDGE *edge, uint64_t numelems) { CBT2(topo, insertEdges, edge, numelems); } int lwt_be_updateEdges(LWT_TOPOLOGY* topo, const LWT_ISO_EDGE* sel_edge, int sel_fields, const LWT_ISO_EDGE* upd_edge, int upd_fields, const LWT_ISO_EDGE* exc_edge, int exc_fields ) { CBT6(topo, updateEdges, sel_edge, sel_fields, upd_edge, upd_fields, exc_edge, exc_fields); } static int lwt_be_updateNodes(LWT_TOPOLOGY* topo, const LWT_ISO_NODE* sel_node, int sel_fields, const LWT_ISO_NODE* upd_node, int upd_fields, const LWT_ISO_NODE* exc_node, int exc_fields ) { CBT6(topo, updateNodes, sel_node, sel_fields, upd_node, upd_fields, exc_node, exc_fields); } static uint64_t lwt_be_updateFacesById(LWT_TOPOLOGY* topo, const LWT_ISO_FACE* faces, uint64_t numfaces ) { CBT2(topo, updateFacesById, faces, numfaces); } static int lwt_be_updateEdgesById(LWT_TOPOLOGY* topo, const LWT_ISO_EDGE* edges, int numedges, int upd_fields ) { CBT3(topo, updateEdgesById, edges, numedges, upd_fields); } static int lwt_be_updateNodesById(LWT_TOPOLOGY* topo, const LWT_ISO_NODE* nodes, int numnodes, int upd_fields ) { CBT3(topo, updateNodesById, nodes, numnodes, upd_fields); } int lwt_be_deleteEdges(LWT_TOPOLOGY* topo, const LWT_ISO_EDGE* sel_edge, int sel_fields ) { CBT2(topo, deleteEdges, sel_edge, sel_fields); } int lwt_be_updateTopoGeomEdgeSplit(LWT_TOPOLOGY* topo, LWT_ELEMID split_edge, LWT_ELEMID new_edge1, LWT_ELEMID new_edge2) { CBT3(topo, updateTopoGeomEdgeSplit, split_edge, new_edge1, new_edge2); } static int lwt_be_updateTopoGeomFaceSplit(LWT_TOPOLOGY* topo, LWT_ELEMID split_face, LWT_ELEMID new_face1, LWT_ELEMID new_face2) { CBT3(topo, updateTopoGeomFaceSplit, split_face, new_face1, new_face2); } static int lwt_be_checkTopoGeomRemEdge(LWT_TOPOLOGY* topo, LWT_ELEMID edge_id, LWT_ELEMID face_left, LWT_ELEMID face_right) { CBT3(topo, checkTopoGeomRemEdge, edge_id, face_left, face_right); } static int lwt_be_checkTopoGeomRemIsoEdge(LWT_TOPOLOGY* topo, LWT_ELEMID edge_id) { CBT1(topo, checkTopoGeomRemIsoEdge, edge_id); } static int lwt_be_checkTopoGeomRemNode(LWT_TOPOLOGY* topo, LWT_ELEMID node_id, LWT_ELEMID eid1, LWT_ELEMID eid2) { CBT3(topo, checkTopoGeomRemNode, node_id, eid1, eid2); } static int lwt_be_checkTopoGeomRemIsoNode(LWT_TOPOLOGY* topo, LWT_ELEMID node_id) { CBT1(topo, checkTopoGeomRemIsoNode, node_id); } static int lwt_be_updateTopoGeomFaceHeal(LWT_TOPOLOGY* topo, LWT_ELEMID face1, LWT_ELEMID face2, LWT_ELEMID newface) { CBT3(topo, updateTopoGeomFaceHeal, face1, face2, newface); } static int lwt_be_updateTopoGeomEdgeHeal(LWT_TOPOLOGY* topo, LWT_ELEMID edge1, LWT_ELEMID edge2, LWT_ELEMID newedge) { CBT3(topo, updateTopoGeomEdgeHeal, edge1, edge2, newedge); } static LWT_ELEMID * lwt_be_getRingEdges(LWT_TOPOLOGY *topo, LWT_ELEMID edge, uint64_t *numedges, uint64_t limit) { CBT3(topo, getRingEdges, edge, numedges, limit); } int lwt_be_ExistsCoincidentNode(LWT_TOPOLOGY* topo, LWPOINT* pt) { uint64_t exists = 0; lwt_be_getNodeWithinDistance2D(topo, pt, 0, &exists, 0, -1); if (exists == UINT64_MAX) { lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return 0; } return exists; } int lwt_be_ExistsEdgeIntersectingPoint(LWT_TOPOLOGY* topo, LWPOINT* pt) { uint64_t exists = 0; lwt_be_getEdgeWithinDistance2D(topo, pt, 0, &exists, 0, -1); if (exists == UINT64_MAX) { lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return 0; } return exists; } static LWT_ISO_EDGE * lwt_be_getClosestEdge(const LWT_TOPOLOGY *topo, const LWPOINT *pt, uint64_t *numelems, int fields) { CBT3(topo, getClosestEdge, pt, numelems, fields); } static GBOX * lwt_be_computeFaceMBR(const LWT_TOPOLOGY *topo, LWT_ELEMID face) { CBT1(topo, computeFaceMBR, face); } /************************************************************************ * * Utility functions * ************************************************************************/ static LWGEOM * _lwt_toposnap(LWGEOM *src, LWGEOM *tgt, double tol) { LWGEOM *tmp = src; LWGEOM *tmp2; int changed; int iterations = 0; int maxiterations = lwgeom_count_vertices(tgt); /* GEOS snapping can be unstable */ /* See https://trac.osgeo.org/geos/ticket/760 */ do { tmp2 = lwgeom_snap(tmp, tgt, tol); ++iterations; changed = ( lwgeom_count_vertices(tmp2) != lwgeom_count_vertices(tmp) ); LWDEBUGF(2, "After iteration %d, geometry changed ? %d (%d vs %d vertices)", iterations, changed, lwgeom_count_vertices(tmp2), lwgeom_count_vertices(tmp)); if ( tmp != src ) lwgeom_free(tmp); tmp = tmp2; } while ( changed && iterations <= maxiterations ); LWDEBUGF(1, "It took %d/%d iterations to properly snap", iterations, maxiterations); return tmp; } static void _lwt_release_faces(LWT_ISO_FACE *faces, int num_faces) { int i; for ( i=0; ibe_iface = iface; topo->be_topo = be_topo; topo->srid = lwt_be_topoGetSRID(topo); topo->hasZ = lwt_be_topoHasZ(topo); topo->precision = lwt_be_topoGetPrecision(topo); return topo; } void lwt_FreeTopology( LWT_TOPOLOGY* topo ) { if ( ! lwt_be_freeTopology(topo) ) { lwnotice("Could not release backend topology memory: %s", lwt_be_lastErrorMessage(topo->be_iface)); } lwfree(topo); } /** * @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; if ( lwpoint_is_empty(pt) ) { lwerror("Cannot add empty point as isolated node"); return -1; } if ( ! skipISOChecks ) { if ( lwt_be_ExistsCoincidentNode(topo, pt) ) /*x*/ { lwerror("SQL/MM Spatial exception - coincident node"); return -1; } if ( lwt_be_ExistsEdgeIntersectingPoint(topo, pt) ) /*x*/ { lwerror("SQL/MM Spatial exception - edge crosses node."); return -1; } } if ( checkFace && ( face == -1 || ! skipISOChecks ) ) { foundInFace = lwt_GetFaceContainingPoint(topo, pt); /*x*/ if ( foundInFace == -1 ) { lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } if ( foundInFace == -1 ) foundInFace = 0; } if ( face == -1 ) { face = foundInFace; } else if ( ! skipISOChecks && foundInFace != face ) { #if 0 lwerror("SQL/MM Spatial exception - within face %d (not %d)", foundInFace, face); #else lwerror("SQL/MM Spatial exception - not within face"); #endif return -1; } LWT_ISO_NODE node; node.node_id = -1; node.containing_face = face; node.geom = pt; if ( ! lwt_be_insertNodes(topo, &node, 1) ) { lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } 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 and * does not have non-boundary intersection with existing edge * * @param myself the id of an edge to skip, if any * (for ChangeEdgeGeom). Can use 0 for none. * * Return -1 on cross or error, 0 if everything is fine. * Note that before returning -1, lwerror is invoked... */ static int _lwt_CheckEdgeCrossing( LWT_TOPOLOGY* topo, LWT_ELEMID start_node, LWT_ELEMID end_node, const LWLINE *geom, LWT_ELEMID myself ) { uint64_t i, num_nodes, num_edges; LWT_ISO_EDGE *edges; LWT_ISO_NODE *nodes; const GBOX *edgebox; GEOSGeometry *edgegg; initGEOS(lwnotice, lwgeom_geos_error); edgegg = LWGEOM2GEOS(lwline_as_lwgeom(geom), 0); if (!edgegg) { lwerror("Could not convert edge geometry to GEOS: %s", lwgeom_geos_errmsg); return -1; } edgebox = lwgeom_get_bbox( lwline_as_lwgeom(geom) ); /* loop over each node within the edge's gbox */ nodes = lwt_be_getNodeWithinBox2D( topo, edgebox, &num_nodes, LWT_COL_NODE_ALL, 0 ); LWDEBUGF(1, "lwt_be_getNodeWithinBox2D returned %d nodes", num_nodes); if (num_nodes == UINT64_MAX) { lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } for ( i=0; inode_id == start_node ) continue; if ( node->node_id == end_node ) continue; /* check if the edge contains this node (not on boundary) */ /* ST_RelateMatch(rec.relate, 'T********') */ pt = getPoint2d_cp(node->geom->point, 0); int contains = ptarray_contains_point_partial(geom->points, pt, LW_FALSE, NULL) == LW_BOUNDARY; if ( contains ) { GEOSGeom_destroy(edgegg); _lwt_release_nodes(nodes, num_nodes); lwerror("SQL/MM Spatial exception - geometry crosses a node"); return -1; } } if ( nodes ) _lwt_release_nodes(nodes, num_nodes); /* may be NULL if num_nodes == 0 */ /* loop over each edge within the edge's gbox */ edges = lwt_be_getEdgeWithinBox2D( topo, edgebox, &num_edges, LWT_COL_EDGE_ALL, 0 ); LWDEBUGF(1, "lwt_be_getEdgeWithinBox2D returned %d edges", num_edges); if (num_edges == UINT64_MAX) { GEOSGeom_destroy(edgegg); lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } for ( i=0; iedge_id; GEOSGeometry *eegg; char *relate; int match; if ( edge_id == myself ) continue; if ( ! edge->geom ) { _lwt_release_edges(edges, num_edges); lwerror("Edge %d has NULL geometry!", edge_id); return -1; } eegg = LWGEOM2GEOS( lwline_as_lwgeom(edge->geom), 0 ); if ( ! eegg ) { GEOSGeom_destroy(edgegg); _lwt_release_edges(edges, num_edges); lwerror("Could not convert edge geometry to GEOS: %s", lwgeom_geos_errmsg); return -1; } LWDEBUGF(2, "Edge %d converted to GEOS", edge_id); /* check if the edge has a non-boundary-boundary intersection with our edge */ relate = GEOSRelateBoundaryNodeRule(eegg, edgegg, 2); if ( ! relate ) { GEOSGeom_destroy(eegg); GEOSGeom_destroy(edgegg); _lwt_release_edges(edges, num_edges); lwerror("GEOSRelateBoundaryNodeRule error: %s", lwgeom_geos_errmsg); return -1; } LWDEBUGF(2, "Edge %d relate pattern is %s", edge_id, relate); match = GEOSRelatePatternMatch(relate, "FF*F*****"); if ( match ) { /* error or no interior intersection */ GEOSGeom_destroy(eegg); GEOSFree(relate); if ( match == 2 ) { _lwt_release_edges(edges, num_edges); GEOSGeom_destroy(edgegg); lwerror("GEOSRelatePatternMatch error: %s", lwgeom_geos_errmsg); return -1; } else continue; /* no interior intersection */ } match = GEOSRelatePatternMatch(relate, "1FFF*FFF2"); if ( match ) { _lwt_release_edges(edges, num_edges); GEOSGeom_destroy(edgegg); GEOSGeom_destroy(eegg); GEOSFree(relate); if ( match == 2 ) { lwerror("GEOSRelatePatternMatch error: %s", lwgeom_geos_errmsg); } else { lwerror("SQL/MM Spatial exception - coincident edge %" LWTFMT_ELEMID, edge_id); } return -1; } match = GEOSRelatePatternMatch(relate, "1********"); if ( match ) { _lwt_release_edges(edges, num_edges); GEOSGeom_destroy(edgegg); GEOSGeom_destroy(eegg); GEOSFree(relate); if ( match == 2 ) { lwerror("GEOSRelatePatternMatch error: %s", lwgeom_geos_errmsg); } else { lwerror("Spatial exception - geometry intersects edge %" LWTFMT_ELEMID, edge_id); } return -1; } match = GEOSRelatePatternMatch(relate, "T********"); if ( match ) { _lwt_release_edges(edges, num_edges); GEOSGeom_destroy(edgegg); GEOSGeom_destroy(eegg); GEOSFree(relate); if ( match == 2 ) { lwerror("GEOSRelatePatternMatch error: %s", lwgeom_geos_errmsg); } else { lwerror("SQL/MM Spatial exception - geometry crosses edge %" LWTFMT_ELEMID, edge_id); } return -1; } match = GEOSRelatePatternMatch(relate, "*T*******"); if ( match ) { _lwt_release_edges(edges, num_edges); GEOSGeom_destroy(edgegg); GEOSGeom_destroy(eegg); GEOSFree(relate); if ( match == 2 ) { lwerror("GEOSRelatePatternMatch error: %s", lwgeom_geos_errmsg); } else { lwerror("Spatial exception - geometry boundary touches interior of edge %" LWTFMT_ELEMID, edge_id); } return -1; } match = GEOSRelatePatternMatch(relate, "***T*****"); if ( match ) { _lwt_release_edges(edges, num_edges); GEOSGeom_destroy(edgegg); GEOSGeom_destroy(eegg); GEOSFree(relate); if ( match == 2 ) { lwerror("GEOSRelatePatternMatch error: %s", lwgeom_geos_errmsg); } else { lwerror("Spatial exception - boundary of edge % touches interior of geometry" LWTFMT_ELEMID, edge_id); } return -1; } LWDEBUGF(2, "Edge %d analisys completed, it does no harm", edge_id); GEOSFree(relate); GEOSGeom_destroy(eegg); } LWDEBUGF(1, "No edge crossing detected among the %d candidate edges", num_edges); if ( edges ) _lwt_release_edges(edges, num_edges); /* would be NULL if num_edges was 0 */ GEOSGeom_destroy(edgegg); return 0; } LWT_ELEMID lwt_AddIsoEdge( LWT_TOPOLOGY* topo, LWT_ELEMID startNode, LWT_ELEMID endNode, const LWLINE* geom ) { uint64_t num_nodes; uint64_t i; LWT_ISO_EDGE newedge; LWT_ISO_NODE *endpoints; LWT_ELEMID containing_face = -1; LWT_ELEMID node_ids[2]; LWT_ISO_NODE updated_nodes[2]; int skipISOChecks = 0; POINT2D p1, p2; /* NOT IN THE SPECS: * A closed edge is never isolated (as it forms a face) */ if (startNode == endNode) { lwerror("Closed edges would not be isolated, try lwt_AddEdgeNewFaces"); return -1; } if ( ! skipISOChecks ) { /* Acurve must be simple */ if ( ! lwgeom_is_simple(lwline_as_lwgeom(geom)) ) { lwerror("SQL/MM Spatial exception - curve not simple"); return -1; } } /* * Check for: * existence of nodes * nodes faces match * Extract: * nodes face id * nodes geoms */ num_nodes = 2; node_ids[0] = startNode; node_ids[1] = endNode; endpoints = lwt_be_getNodeById( topo, node_ids, &num_nodes, LWT_COL_NODE_ALL ); if (num_nodes == UINT64_MAX) { lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } else if ( num_nodes < 2 ) { if ( num_nodes ) _lwt_release_nodes(endpoints, num_nodes); lwerror("SQL/MM Spatial exception - non-existent node"); return -1; } for ( i=0; icontaining_face == -1 ) { _lwt_release_nodes(endpoints, num_nodes); lwerror("SQL/MM Spatial exception - not isolated node"); return -1; } if ( containing_face == -1 ) containing_face = n->containing_face; else if ( containing_face != n->containing_face ) { _lwt_release_nodes(endpoints, num_nodes); lwerror("SQL/MM Spatial exception - nodes in different faces"); return -1; } if ( ! skipISOChecks ) { if ( n->node_id == startNode ) { /* l) Check that start point of acurve match start node geoms. */ getPoint2d_p(geom->points, 0, &p1); getPoint2d_p(n->geom->point, 0, &p2); if ( ! P2D_SAME_STRICT(&p1, &p2) ) { _lwt_release_nodes(endpoints, num_nodes); lwerror("SQL/MM Spatial exception - " "start node not geometry start point."); return -1; } } else { /* m) Check that end point of acurve match end node geoms. */ getPoint2d_p(geom->points, geom->points->npoints-1, &p1); getPoint2d_p(n->geom->point, 0, &p2); if ( ! P2D_SAME_STRICT(&p1, &p2) ) { _lwt_release_nodes(endpoints, num_nodes); lwerror("SQL/MM Spatial exception - " "end node not geometry end point."); return -1; } } } } if ( num_nodes ) _lwt_release_nodes(endpoints, num_nodes); if ( ! skipISOChecks ) { if ( _lwt_CheckEdgeCrossing( topo, startNode, endNode, geom, 0 ) ) { /* would have called lwerror already, leaking :( */ return -1; } } /* * All checks passed, time to prepare the new edge */ newedge.edge_id = lwt_be_getNextEdgeId( topo ); if ( newedge.edge_id == -1 ) { lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } /* TODO: this should likely be an exception instead ! */ if ( containing_face == -1 ) containing_face = 0; newedge.start_node = startNode; newedge.end_node = endNode; newedge.face_left = newedge.face_right = containing_face; newedge.next_left = -newedge.edge_id; newedge.next_right = newedge.edge_id; newedge.geom = (LWLINE *)geom; /* const cast.. */ int ret = lwt_be_insertEdges(topo, &newedge, 1); if ( ret == -1 ) { lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } else if ( ret == 0 ) { lwerror("Insertion of split edge failed (no reason)"); return -1; } /* * Update Node containing_face values * * the nodes anode and anothernode are no more isolated * because now there is an edge connecting them */ updated_nodes[0].node_id = startNode; updated_nodes[0].containing_face = -1; updated_nodes[1].node_id = endNode; updated_nodes[1].containing_face = -1; ret = lwt_be_updateNodesById(topo, updated_nodes, 2, LWT_COL_NODE_CONTAINING_FACE); if ( ret == -1 ) { lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } return newedge.edge_id; } static LWCOLLECTION * _lwt_EdgeSplit( LWT_TOPOLOGY* topo, LWT_ELEMID edge, LWPOINT* pt, int skipISOChecks, LWT_ISO_EDGE** oldedge ) { LWGEOM *split; LWCOLLECTION *split_col; uint64_t i; /* Get edge */ i = 1; LWDEBUG(1, "calling lwt_be_getEdgeById"); *oldedge = lwt_be_getEdgeById(topo, &edge, &i, LWT_COL_EDGE_ALL); LWDEBUGF(1, "lwt_be_getEdgeById returned %p", *oldedge); if ( ! *oldedge ) { LWDEBUGF(1, "lwt_be_getEdgeById returned NULL and set i=%d", i); if (i == UINT64_MAX) { lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return NULL; } else if ( i == 0 ) { lwerror("SQL/MM Spatial exception - non-existent edge"); return NULL; } else { lwerror("Backend coding error: getEdgeById callback returned NULL " "but numelements output parameter has value %d " "(expected 0 or 1)", i); return NULL; } } /* * - check if a coincident node already exists */ if ( ! skipISOChecks ) { LWDEBUG(1, "calling lwt_be_ExistsCoincidentNode"); if ( lwt_be_ExistsCoincidentNode(topo, pt) ) /*x*/ { LWDEBUG(1, "lwt_be_ExistsCoincidentNode returned"); _lwt_release_edges(*oldedge, 1); lwerror("SQL/MM Spatial exception - coincident node"); return NULL; } LWDEBUG(1, "lwt_be_ExistsCoincidentNode returned"); } /* Split edge */ split = lwgeom_split((LWGEOM*)(*oldedge)->geom, (LWGEOM*)pt); if ( ! split ) { _lwt_release_edges(*oldedge, 1); lwerror("could not split edge by point ?"); return NULL; } split_col = lwgeom_as_lwcollection(split); if ( ! split_col ) { _lwt_release_edges(*oldedge, 1); lwgeom_free(split); lwerror("lwgeom_as_lwcollection returned NULL"); return NULL; } if (split_col->ngeoms < 2) { _lwt_release_edges(*oldedge, 1); lwgeom_free(split); lwerror("SQL/MM Spatial exception - point not on edge"); return NULL; } #if 0 { size_t sz; char *wkt = lwgeom_to_wkt((LWGEOM*)split_col, WKT_EXTENDED, 2, &sz); LWDEBUGF(1, "returning split col: %s", wkt); lwfree(wkt); } #endif return split_col; } LWT_ELEMID lwt_ModEdgeSplit( LWT_TOPOLOGY* topo, LWT_ELEMID edge, LWPOINT* pt, int skipISOChecks ) { LWT_ISO_NODE node; LWT_ISO_EDGE* oldedge = NULL; LWCOLLECTION *split_col; const LWGEOM *oldedge_geom; const LWGEOM *newedge_geom; LWT_ISO_EDGE newedge1; LWT_ISO_EDGE seledge, updedge, excedge; int ret; split_col = _lwt_EdgeSplit( topo, edge, pt, skipISOChecks, &oldedge ); if ( ! split_col ) return -1; /* should have raised an exception */ oldedge_geom = split_col->geoms[0]; newedge_geom = split_col->geoms[1]; /* Make sure the SRID is set on the subgeom */ ((LWGEOM*)oldedge_geom)->srid = split_col->srid; ((LWGEOM*)newedge_geom)->srid = split_col->srid; /* Add new node, getting new id back */ node.node_id = -1; node.containing_face = -1; /* means not-isolated */ node.geom = pt; if ( ! lwt_be_insertNodes(topo, &node, 1) ) { _lwt_release_edges(oldedge, 1); lwcollection_free(split_col); lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } if (node.node_id == -1) { /* should have been set by backend */ _lwt_release_edges(oldedge, 1); lwcollection_free(split_col); lwerror("Backend coding error: " "insertNodes callback did not return node_id"); return -1; } /* Insert the new edge */ newedge1.edge_id = lwt_be_getNextEdgeId(topo); if ( newedge1.edge_id == -1 ) { _lwt_release_edges(oldedge, 1); lwcollection_free(split_col); lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } newedge1.start_node = node.node_id; newedge1.end_node = oldedge->end_node; newedge1.face_left = oldedge->face_left; newedge1.face_right = oldedge->face_right; newedge1.next_left = oldedge->next_left == -oldedge->edge_id ? -newedge1.edge_id : oldedge->next_left; newedge1.next_right = -oldedge->edge_id; newedge1.geom = lwgeom_as_lwline(newedge_geom); /* lwgeom_split of a line should only return lines ... */ if ( ! newedge1.geom ) { _lwt_release_edges(oldedge, 1); lwcollection_free(split_col); lwerror("first geometry in lwgeom_split output is not a line"); return -1; } ret = lwt_be_insertEdges(topo, &newedge1, 1); if ( ret == -1 ) { _lwt_release_edges(oldedge, 1); lwcollection_free(split_col); lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } else if ( ret == 0 ) { _lwt_release_edges(oldedge, 1); lwcollection_free(split_col); lwerror("Insertion of split edge failed (no reason)"); return -1; } /* Update the old edge */ updedge.geom = lwgeom_as_lwline(oldedge_geom); /* lwgeom_split of a line should only return lines ... */ if ( ! updedge.geom ) { _lwt_release_edges(oldedge, 1); lwcollection_free(split_col); lwerror("second geometry in lwgeom_split output is not a line"); return -1; } updedge.next_left = newedge1.edge_id; updedge.end_node = node.node_id; ret = lwt_be_updateEdges(topo, oldedge, LWT_COL_EDGE_EDGE_ID, &updedge, LWT_COL_EDGE_GEOM|LWT_COL_EDGE_NEXT_LEFT|LWT_COL_EDGE_END_NODE, NULL, 0); if ( ret == -1 ) { _lwt_release_edges(oldedge, 1); lwcollection_free(split_col); lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } else if ( ret == 0 ) { _lwt_release_edges(oldedge, 1); lwcollection_free(split_col); lwerror("Edge being split (%d) disappeared during operations?", oldedge->edge_id); return -1; } else if ( ret > 1 ) { _lwt_release_edges(oldedge, 1); lwcollection_free(split_col); lwerror("More than a single edge found with id %d !", oldedge->edge_id); return -1; } /* Update all next edge references to match new layout (ST_ModEdgeSplit) */ updedge.next_right = -newedge1.edge_id; excedge.edge_id = newedge1.edge_id; seledge.next_right = -oldedge->edge_id; seledge.start_node = oldedge->end_node; ret = lwt_be_updateEdges(topo, &seledge, LWT_COL_EDGE_NEXT_RIGHT|LWT_COL_EDGE_START_NODE, &updedge, LWT_COL_EDGE_NEXT_RIGHT, &excedge, LWT_COL_EDGE_EDGE_ID); if ( ret == -1 ) { _lwt_release_edges(oldedge, 1); lwcollection_free(split_col); lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } updedge.next_left = -newedge1.edge_id; excedge.edge_id = newedge1.edge_id; seledge.next_left = -oldedge->edge_id; seledge.end_node = oldedge->end_node; ret = lwt_be_updateEdges(topo, &seledge, LWT_COL_EDGE_NEXT_LEFT|LWT_COL_EDGE_END_NODE, &updedge, LWT_COL_EDGE_NEXT_LEFT, &excedge, LWT_COL_EDGE_EDGE_ID); if ( ret == -1 ) { _lwt_release_edges(oldedge, 1); lwcollection_free(split_col); lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } /* Update TopoGeometries composition */ ret = lwt_be_updateTopoGeomEdgeSplit(topo, oldedge->edge_id, newedge1.edge_id, -1); if ( ! ret ) { _lwt_release_edges(oldedge, 1); lwcollection_free(split_col); lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } _lwt_release_edges(oldedge, 1); lwcollection_free(split_col); /* return new node id */ return node.node_id; } LWT_ELEMID lwt_NewEdgesSplit( LWT_TOPOLOGY* topo, LWT_ELEMID edge, LWPOINT* pt, int skipISOChecks ) { LWT_ISO_NODE node; LWT_ISO_EDGE* oldedge = NULL; LWCOLLECTION *split_col; const LWGEOM *oldedge_geom; const LWGEOM *newedge_geom; LWT_ISO_EDGE newedges[2]; LWT_ISO_EDGE seledge, updedge; int ret; split_col = _lwt_EdgeSplit( topo, edge, pt, skipISOChecks, &oldedge ); if ( ! split_col ) return -1; /* should have raised an exception */ oldedge_geom = split_col->geoms[0]; newedge_geom = split_col->geoms[1]; /* Make sure the SRID is set on the subgeom */ ((LWGEOM*)oldedge_geom)->srid = split_col->srid; ((LWGEOM*)newedge_geom)->srid = split_col->srid; /* Add new node, getting new id back */ node.node_id = -1; node.containing_face = -1; /* means not-isolated */ node.geom = pt; if ( ! lwt_be_insertNodes(topo, &node, 1) ) { _lwt_release_edges(oldedge, 1); lwcollection_free(split_col); lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } if (node.node_id == -1) { _lwt_release_edges(oldedge, 1); lwcollection_free(split_col); /* should have been set by backend */ lwerror("Backend coding error: " "insertNodes callback did not return node_id"); return -1; } /* Delete the old edge */ seledge.edge_id = edge; ret = lwt_be_deleteEdges(topo, &seledge, LWT_COL_EDGE_EDGE_ID); if ( ret == -1 ) { _lwt_release_edges(oldedge, 1); lwcollection_free(split_col); lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } /* Get new edges identifiers */ newedges[0].edge_id = lwt_be_getNextEdgeId(topo); if ( newedges[0].edge_id == -1 ) { _lwt_release_edges(oldedge, 1); lwcollection_free(split_col); lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } newedges[1].edge_id = lwt_be_getNextEdgeId(topo); if ( newedges[1].edge_id == -1 ) { _lwt_release_edges(oldedge, 1); lwcollection_free(split_col); lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } /* Define the first new edge (to new node) */ newedges[0].start_node = oldedge->start_node; newedges[0].end_node = node.node_id; newedges[0].face_left = oldedge->face_left; newedges[0].face_right = oldedge->face_right; newedges[0].next_left = newedges[1].edge_id; if ( oldedge->next_right == edge ) newedges[0].next_right = newedges[0].edge_id; else if ( oldedge->next_right == -edge ) newedges[0].next_right = -newedges[1].edge_id; else newedges[0].next_right = oldedge->next_right; newedges[0].geom = lwgeom_as_lwline(oldedge_geom); /* lwgeom_split of a line should only return lines ... */ if ( ! newedges[0].geom ) { _lwt_release_edges(oldedge, 1); lwcollection_free(split_col); lwerror("first geometry in lwgeom_split output is not a line"); return -1; } /* Define the second new edge (from new node) */ newedges[1].start_node = node.node_id; newedges[1].end_node = oldedge->end_node; newedges[1].face_left = oldedge->face_left; newedges[1].face_right = oldedge->face_right; newedges[1].next_right = -newedges[0].edge_id; if ( oldedge->next_left == -edge ) newedges[1].next_left = -newedges[1].edge_id; else if ( oldedge->next_left == edge ) newedges[1].next_left = newedges[0].edge_id; else newedges[1].next_left = oldedge->next_left; newedges[1].geom = lwgeom_as_lwline(newedge_geom); /* lwgeom_split of a line should only return lines ... */ if ( ! newedges[1].geom ) { _lwt_release_edges(oldedge, 1); lwcollection_free(split_col); lwerror("second geometry in lwgeom_split output is not a line"); return -1; } /* Insert both new edges */ ret = lwt_be_insertEdges(topo, newedges, 2); if ( ret == -1 ) { _lwt_release_edges(oldedge, 1); lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } else if ( ret == 0 ) { _lwt_release_edges(oldedge, 1); lwcollection_free(split_col); lwerror("Insertion of split edge failed (no reason)"); return -1; } /* Update all next edge references pointing to old edge id */ updedge.next_right = newedges[1].edge_id; seledge.next_right = edge; seledge.start_node = oldedge->start_node; ret = lwt_be_updateEdges(topo, &seledge, LWT_COL_EDGE_NEXT_RIGHT|LWT_COL_EDGE_START_NODE, &updedge, LWT_COL_EDGE_NEXT_RIGHT, NULL, 0); if ( ret == -1 ) { _lwt_release_edges(oldedge, 1); lwcollection_free(split_col); lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } updedge.next_right = -newedges[0].edge_id; seledge.next_right = -edge; seledge.start_node = oldedge->end_node; ret = lwt_be_updateEdges(topo, &seledge, LWT_COL_EDGE_NEXT_RIGHT|LWT_COL_EDGE_START_NODE, &updedge, LWT_COL_EDGE_NEXT_RIGHT, NULL, 0); if ( ret == -1 ) { _lwt_release_edges(oldedge, 1); lwcollection_free(split_col); lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } updedge.next_left = newedges[0].edge_id; seledge.next_left = edge; seledge.end_node = oldedge->start_node; ret = lwt_be_updateEdges(topo, &seledge, LWT_COL_EDGE_NEXT_LEFT|LWT_COL_EDGE_END_NODE, &updedge, LWT_COL_EDGE_NEXT_LEFT, NULL, 0); if ( ret == -1 ) { _lwt_release_edges(oldedge, 1); lwcollection_free(split_col); lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } updedge.next_left = -newedges[1].edge_id; seledge.next_left = -edge; seledge.end_node = oldedge->end_node; ret = lwt_be_updateEdges(topo, &seledge, LWT_COL_EDGE_NEXT_LEFT|LWT_COL_EDGE_END_NODE, &updedge, LWT_COL_EDGE_NEXT_LEFT, NULL, 0); if ( ret == -1 ) { _lwt_release_edges(oldedge, 1); lwcollection_release(split_col); lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } /* Update TopoGeometries composition */ ret = lwt_be_updateTopoGeomEdgeSplit(topo, oldedge->edge_id, newedges[0].edge_id, newedges[1].edge_id); if ( ! ret ) { _lwt_release_edges(oldedge, 1); lwcollection_free(split_col); lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } _lwt_release_edges(oldedge, 1); lwcollection_free(split_col); /* return new node id */ return node.node_id; } /* Data structure used by AddEdgeX functions */ typedef struct edgeend_t { /* Signed identifier of next clockwise edge (+outgoing,-incoming) */ LWT_ELEMID nextCW; /* Identifier of face between myaz and next CW edge */ LWT_ELEMID cwFace; /* Signed identifier of next counterclockwise edge (+outgoing,-incoming) */ LWT_ELEMID nextCCW; /* Identifier of face between myaz and next CCW edge */ LWT_ELEMID ccwFace; int was_isolated; double myaz; /* azimuth of edgeend geometry */ } edgeend; /* * Get first distinct vertex from endpoint * @param pa the pointarray to seek points in * @param ref the point we want to search a distinct one * @param from vertex index to start from (will really start from "from"+dir) * @param dir 1 to go forward * -1 to go backward * @return 0 if edge is collapsed (no distinct points) */ static int _lwt_FirstDistinctVertex2D(const POINTARRAY* pa, POINT2D *ref, int from, int dir, POINT2D *op) { int i, toofar, inc; POINT2D fp; if ( dir > 0 ) { toofar = pa->npoints; inc = 1; } else { toofar = -1; inc = -1; } LWDEBUGF(1, "first point is index %d", from); fp = *ref; /* getPoint2d_p(pa, from, &fp); */ for ( i = from+inc; i != toofar; i += inc ) { LWDEBUGF(1, "testing point %d", i); getPoint2d_p(pa, i, op); /* pick next point */ if ( P2D_SAME_STRICT(op,&fp) ) continue; /* equal to startpoint */ /* this is a good one, neither same of start nor of end point */ return 1; /* found */ } /* no distinct vertices found */ return 0; } /* * Return non-zero on failure (lwerror is invoked in that case) * Possible failures: * -1 no two distinct vertices exist * -2 azimuth computation failed for first edge end */ static int _lwt_InitEdgeEndByLine(edgeend *fee, edgeend *lee, LWLINE *edge, POINT2D *fp, POINT2D *lp) { POINTARRAY *pa = edge->points; POINT2D pt; fee->nextCW = fee->nextCCW = lee->nextCW = lee->nextCCW = 0; fee->cwFace = fee->ccwFace = lee->cwFace = lee->ccwFace = -1; /* Compute azimuth of first edge end */ LWDEBUG(1, "computing azimuth of first edge end"); if ( ! _lwt_FirstDistinctVertex2D(pa, fp, 0, 1, &pt) ) { lwerror("Invalid edge (no two distinct vertices exist)"); return -1; } if ( ! azimuth_pt_pt(fp, &pt, &(fee->myaz)) ) { lwerror("error computing azimuth of first edgeend [%.15g %.15g,%.15g %.15g]", fp->x, fp->y, pt.x, pt.y); return -2; } LWDEBUGF(1, "azimuth of first edge end [%.15g %.15g,%.15g %.15g] is %g", fp->x, fp->y, pt.x, pt.y, fee->myaz); /* Compute azimuth of second edge end */ LWDEBUG(1, "computing azimuth of second edge end"); if ( ! _lwt_FirstDistinctVertex2D(pa, lp, pa->npoints-1, -1, &pt) ) { lwerror("Invalid edge (no two distinct vertices exist)"); return -1; } if ( ! azimuth_pt_pt(lp, &pt, &(lee->myaz)) ) { lwerror("error computing azimuth of last edgeend [%.15g %.15g,%.15g %.15g]", lp->x, lp->y, pt.x, pt.y); return -2; } LWDEBUGF(1, "azimuth of last edge end [%.15g %.15g,%.15g %.15g] is %g", lp->x, lp->y, pt.x, pt.y, lee->myaz); return 0; } /* * Find the first edges encountered going clockwise and counterclockwise * around a node, starting from the given azimuth, and take * note of the face on the both sides. * * @param topo the topology to act upon * @param node the identifier of the node to analyze * @param data input (myaz) / output (nextCW, nextCCW) parameter * @param other edgeend, if also incident to given node (closed edge). * @param myedge_id identifier of the edge id that data->myaz belongs to * @return number of incident edges found * */ static int _lwt_FindAdjacentEdges( LWT_TOPOLOGY* topo, LWT_ELEMID node, edgeend *data, edgeend *other, int myedge_id ) { LWT_ISO_EDGE *edges; uint64_t numedges = 1; uint64_t i; double minaz, maxaz; double az, azdif; data->nextCW = data->nextCCW = 0; data->cwFace = data->ccwFace = -1; if ( other ) { azdif = other->myaz - data->myaz; if ( azdif < 0 ) azdif += 2 * M_PI; minaz = maxaz = azdif; /* TODO: set nextCW/nextCCW/cwFace/ccwFace to other->something ? */ LWDEBUGF(1, "Other edge end has cwFace=%d and ccwFace=%d", other->cwFace, other->ccwFace); } else { minaz = maxaz = -1; } LWDEBUGF(1, "Looking for edges incident to node %" LWTFMT_ELEMID " and adjacent to azimuth %g", node, data->myaz); /* Get incident edges */ edges = lwt_be_getEdgeByNode( topo, &node, &numedges, LWT_COL_EDGE_ALL ); if (numedges == UINT64_MAX) { lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return 0; } LWDEBUGF(1, "getEdgeByNode returned %d edges, minaz=%g, maxaz=%g", numedges, minaz, maxaz); /* For each incident edge-end (1 or 2): */ for ( i = 0; i < numedges; ++i ) { LWT_ISO_EDGE *edge; LWGEOM *g; LWGEOM *cleangeom; POINT2D p1, p2; POINTARRAY *pa; edge = &(edges[i]); if ( edge->edge_id == myedge_id ) continue; g = lwline_as_lwgeom(edge->geom); /* NOTE: remove_repeated_points call could be replaced by * some other mean to pick two distinct points for endpoints */ cleangeom = lwgeom_remove_repeated_points( g, 0 ); pa = lwgeom_as_lwline(cleangeom)->points; if ( pa->npoints < 2 ) {{ LWT_ELEMID id = edge->edge_id; _lwt_release_edges(edges, numedges); lwgeom_free(cleangeom); lwerror("corrupted topology: edge %" LWTFMT_ELEMID " does not have two distinct points", id); return -1; }} if ( edge->start_node == node ) { getPoint2d_p(pa, 0, &p1); if ( ! _lwt_FirstDistinctVertex2D(pa, &p1, 0, 1, &p2) ) { lwerror("Edge %d has no distinct vertices: [%.15g %.15g,%.15g %.15g]: ", edge->edge_id, p1.x, p1.y, p2.x, p2.y); return -1; } LWDEBUGF(1, "edge %" LWTFMT_ELEMID " starts on node %" LWTFMT_ELEMID ", edgeend is [%.15g %.15g,%.15g %.15g]", edge->edge_id, node, p1.x, p1.y, p2.x, p2.y); if ( ! azimuth_pt_pt(&p1, &p2, &az) ) {{ LWT_ELEMID id = edge->edge_id; _lwt_release_edges(edges, numedges); lwgeom_free(cleangeom); lwerror("error computing azimuth of edge %d first edgeend [%.15g %.15g,%.15g %.15g]", id, p1.x, p1.y, p2.x, p2.y); return -1; }} azdif = az - data->myaz; LWDEBUGF(1, "azimuth of edge %" LWTFMT_ELEMID ": %.15g (diff: %.15g)", edge->edge_id, az, azdif); if ( azdif < 0 ) azdif += 2 * M_PI; if ( minaz == -1 ) { minaz = maxaz = azdif; data->nextCW = data->nextCCW = edge->edge_id; /* outgoing */ data->cwFace = edge->face_left; data->ccwFace = edge->face_right; LWDEBUGF(1, "new nextCW and nextCCW edge is %" LWTFMT_ELEMID ", outgoing, " "with face_left %" LWTFMT_ELEMID " and face_right %" LWTFMT_ELEMID " (face_right is new ccwFace, face_left is new cwFace)", edge->edge_id, edge->face_left, edge->face_right); } else { if ( azdif < minaz ) { data->nextCW = edge->edge_id; /* outgoing */ data->cwFace = edge->face_left; LWDEBUGF(1, "new nextCW edge is %" LWTFMT_ELEMID ", outgoing, " "with face_left %" LWTFMT_ELEMID " and face_right %" LWTFMT_ELEMID " (previous had minaz=%g, face_left is new cwFace)", edge->edge_id, edge->face_left, edge->face_right, minaz); minaz = azdif; } else if ( azdif > maxaz ) { data->nextCCW = edge->edge_id; /* outgoing */ data->ccwFace = edge->face_right; LWDEBUGF(1, "new nextCCW edge is %" LWTFMT_ELEMID ", outgoing, " "with face_left %" LWTFMT_ELEMID " and face_right %" LWTFMT_ELEMID " (previous had maxaz=%g, face_right is new ccwFace)", edge->edge_id, edge->face_left, edge->face_right, maxaz); maxaz = azdif; } } } if ( edge->end_node == node ) { getPoint2d_p(pa, pa->npoints-1, &p1); if ( ! _lwt_FirstDistinctVertex2D(pa, &p1, pa->npoints-1, -1, &p2) ) { lwerror("Edge %d has no distinct vertices: [%.15g %.15g,%.15g %.15g]: ", edge->edge_id, p1.x, p1.y, p2.x, p2.y); return -1; } LWDEBUGF(1, "edge %" LWTFMT_ELEMID " ends on node %" LWTFMT_ELEMID ", edgeend is [%.15g %.15g,%.15g %.15g]", edge->edge_id, node, p1.x, p1.y, p2.x, p2.y); if ( ! azimuth_pt_pt(&p1, &p2, &az) ) {{ LWT_ELEMID id = edge->edge_id; _lwt_release_edges(edges, numedges); lwgeom_free(cleangeom); lwerror("error computing azimuth of edge %d last edgeend [%.15g %.15g,%.15g %.15g]", id, p1.x, p1.y, p2.x, p2.y); return -1; }} azdif = az - data->myaz; LWDEBUGF(1, "azimuth of edge %" LWTFMT_ELEMID ": %.15g (diff: %.15g)", edge->edge_id, az, azdif); if ( azdif < 0 ) azdif += 2 * M_PI; if ( minaz == -1 ) { minaz = maxaz = azdif; data->nextCW = data->nextCCW = -edge->edge_id; /* incoming */ data->cwFace = edge->face_right; data->ccwFace = edge->face_left; LWDEBUGF(1, "new nextCW and nextCCW edge is %" LWTFMT_ELEMID ", incoming, " "with face_left %" LWTFMT_ELEMID " and face_right %" LWTFMT_ELEMID " (face_right is new cwFace, face_left is new ccwFace)", edge->edge_id, edge->face_left, edge->face_right); } else { if ( azdif < minaz ) { data->nextCW = -edge->edge_id; /* incoming */ data->cwFace = edge->face_right; LWDEBUGF(1, "new nextCW edge is %" LWTFMT_ELEMID ", incoming, " "with face_left %" LWTFMT_ELEMID " and face_right %" LWTFMT_ELEMID " (previous had minaz=%g, face_right is new cwFace)", edge->edge_id, edge->face_left, edge->face_right, minaz); minaz = azdif; } else if ( azdif > maxaz ) { data->nextCCW = -edge->edge_id; /* incoming */ data->ccwFace = edge->face_left; LWDEBUGF(1, "new nextCCW edge is %" LWTFMT_ELEMID ", outgoing, from start point, " "with face_left %" LWTFMT_ELEMID " and face_right %" LWTFMT_ELEMID " (previous had maxaz=%g, face_left is new ccwFace)", edge->edge_id, edge->face_left, edge->face_right, maxaz); maxaz = azdif; } } } lwgeom_free(cleangeom); } if ( numedges ) _lwt_release_edges(edges, numedges); LWDEBUGF(1, "edges adjacent to azimuth %g" " (incident to node %" LWTFMT_ELEMID ")" ": CW:%" LWTFMT_ELEMID "(%g) CCW:%" LWTFMT_ELEMID "(%g)", data->myaz, node, data->nextCW, minaz, data->nextCCW, maxaz); if ( myedge_id < 1 && numedges && data->cwFace != data->ccwFace ) { if ( data->cwFace != -1 && data->ccwFace != -1 ) { lwerror("Corrupted topology: adjacent edges %" LWTFMT_ELEMID " and %" LWTFMT_ELEMID " bind different face (%" LWTFMT_ELEMID " and %" LWTFMT_ELEMID ")", data->nextCW, data->nextCCW, data->cwFace, data->ccwFace); return -1; } } /* Return number of incident edges found */ return numedges; } /* * Get a point internal to the line and write it into the "ip" * parameter * * return 0 on failure (line is empty or collapsed), 1 otherwise */ static int _lwt_GetInteriorEdgePoint(const LWLINE* edge, POINT2D* ip) { uint32_t i; POINT2D fp, lp, tp; POINTARRAY *pa = edge->points; if ( pa->npoints < 2 ) return 0; /* empty or structurally collapsed */ getPoint2d_p(pa, 0, &fp); /* save first point */ getPoint2d_p(pa, pa->npoints-1, &lp); /* save last point */ for (i=1; inpoints-1; ++i) { getPoint2d_p(pa, i, &tp); /* pick next point */ if ( P2D_SAME_STRICT(&tp, &fp) ) continue; /* equal to startpoint */ if ( P2D_SAME_STRICT(&tp, &lp) ) continue; /* equal to endpoint */ /* this is a good one, neither same of start nor of end point */ *ip = tp; return 1; /* found */ } /* no distinct vertex found */ /* interpolate if start point != end point */ if ( P2D_SAME_STRICT(&fp, &lp) ) return 0; /* no distinct points in edge */ ip->x = fp.x + ( (lp.x - fp.x) * 0.5 ); ip->y = fp.y + ( (lp.y - fp.y) * 0.5 ); return 1; } static LWPOLY * _lwt_MakeRingShell(LWT_TOPOLOGY *topo, LWT_ELEMID *signed_edge_ids, uint64_t num_signed_edge_ids) { LWT_ELEMID *edge_ids; uint64_t numedges, i, j; LWT_ISO_EDGE *ring_edges; /* Construct a polygon using edges of the ring */ numedges = 0; edge_ids = lwalloc(sizeof(LWT_ELEMID)*num_signed_edge_ids); for (i=0; ibe_iface)); return NULL; } else if ( i != numedges ) { lwfree( signed_edge_ids ); _lwt_release_edges(ring_edges, i); lwerror("Unexpected error: %d edges found when expecting %d", i, numedges); return NULL; } /* Should now build a polygon with those edges, in the order * given by GetRingEdges. */ POINTARRAY *pa = NULL; for ( i=0; igeom->points); if ( eid < 0 ) ptarray_reverse_in_place(pa); } else { if ( eid < 0 ) { epa = ptarray_clone_deep(edge->geom->points); ptarray_reverse_in_place(epa); ptarray_append_ptarray(pa, epa, 0); ptarray_free(epa); } else { /* avoid a clone here */ ptarray_append_ptarray(pa, edge->geom->points, 0); } } } _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 ) { uint64_t numfaceedges, i, j; int newface_outside; uint64_t 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]; if ( ! ptarray_is_closed_2d(pa) ) { lwpoly_free(shell); lwfree( signed_edge_ids ); lwerror("Corrupted topology: ring of edge %" LWTFMT_ELEMID " is geometrically not-closed", sedge); return -2; } int isccw = ptarray_isccw(pa); LWDEBUGF(1, "Ring of edge %" LWTFMT_ELEMID " is %sclockwise", sedge, isccw ? "counter" : ""); const GBOX* shellbox = lwgeom_get_bbox(lwpoly_as_lwgeom(shell)); if ( face == 0 ) { /* Edge split the universe face */ if ( ! isccw ) { lwpoly_free(shell); lwfree( signed_edge_ids ); /* Face on the left side of this ring is the universe face. * Next call (for the other side) should create the split face */ LWDEBUG(1, "The left face of this clockwise ring is the universe, " "won't create a new face there"); return -1; } } if ( mbr_only && face != 0 ) { if ( isccw ) {{ LWT_ISO_FACE updface; updface.face_id = face; updface.mbr = (GBOX *)shellbox; /* const cast, we won't free it, later */ int ret = lwt_be_updateFacesById( topo, &updface, 1 ); if ( ret == -1 ) { lwfree( signed_edge_ids ); lwpoly_free(shell); /* NOTE: owns shellbox above */ lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -2; } if ( ret != 1 ) { lwfree( signed_edge_ids ); lwpoly_free(shell); /* NOTE: owns shellbox above */ lwerror("Unexpected error: %d faces found when expecting 1", ret); return -2; } }} lwfree( signed_edge_ids ); lwpoly_free(shell); /* NOTE: owns shellbox above */ return -1; /* mbr only was requested */ } LWT_ISO_FACE *oldface = NULL; LWT_ISO_FACE newface; newface.face_id = -1; if ( face != 0 && ! isccw) {{ /* Face created an hole in an outer face */ uint64_t nfaces = 1; oldface = lwt_be_getFaceById(topo, &face, &nfaces, LWT_COL_FACE_ALL); if (nfaces == UINT64_MAX) { lwfree( signed_edge_ids ); lwpoly_free(shell); /* NOTE: owns shellbox */ lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -2; } if ( nfaces != 1 ) { lwfree( signed_edge_ids ); lwpoly_free(shell); /* NOTE: owns shellbox */ lwerror("Unexpected error: %d faces found when expecting 1", nfaces); return -2; } newface.mbr = oldface->mbr; }} else { newface.mbr = (GBOX *)shellbox; /* const cast, we won't free it, later */ } /* Insert the new face */ int ret = lwt_be_insertFaces( topo, &newface, 1 ); if ( ret == -1 ) { lwfree( signed_edge_ids ); lwpoly_free(shell); /* NOTE: owns shellbox */ lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -2; } if ( ret != 1 ) { lwfree( signed_edge_ids ); lwpoly_free(shell); /* NOTE: owns shellbox */ lwerror("Unexpected error: %d faces inserted when expecting 1", ret); return -2; } if ( oldface ) { newface.mbr = NULL; /* it is a reference to oldface mbr... */ _lwt_release_faces(oldface, 1); } /* Update side location of new face edges */ /* We want the new face to be on the left, if possible */ if ( face != 0 && ! isccw ) { /* ring is clockwise in a real face */ /* face shrinked, must update all non-contained edges and nodes */ LWDEBUG(1, "New face is on the outside of the ring, updating rings in former shell"); newface_outside = 1; /* newface is outside */ } else { LWDEBUG(1, "New face is on the inside of the ring, updating forward edges in new ring"); newface_outside = 0; /* newface is inside */ } /* Update edges bounding the old face */ /* (1) fetch all edges where left_face or right_face is = oldface */ int fields = LWT_COL_EDGE_EDGE_ID | LWT_COL_EDGE_FACE_LEFT | LWT_COL_EDGE_FACE_RIGHT | LWT_COL_EDGE_GEOM ; numfaceedges = 1; edges = lwt_be_getEdgeByFace( topo, &face, &numfaceedges, fields, newface.mbr ); if (numfaceedges == UINT64_MAX) { lwfree(signed_edge_ids); lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -2; } LWDEBUGF(1, "_lwt_AddFaceSplit: lwt_be_getEdgeByFace(%d) returned %d edges", face, numfaceedges); if ( numfaceedges ) { forward_edges = lwalloc(sizeof(LWT_ISO_EDGE)*numfaceedges); forward_edges_count = 0; backward_edges = lwalloc(sizeof(LWT_ISO_EDGE)*numfaceedges); backward_edges_count = 0; /* (2) loop over the results and: */ for ( i=0; iedge_id ) { /* IDEA: remove entry from signed_edge_ids, to speed next loop ? */ LWDEBUGF(1, "Edge %d is a known forward edge of the new ring", e->edge_id); forward_edges[forward_edges_count].edge_id = e->edge_id; forward_edges[forward_edges_count++].face_left = newface.face_id; found++; if ( found == 2 ) break; /* both edge sides are found on the ring */ } else if ( -seid == e->edge_id ) { /* IDEA: remove entry from signed_edge_ids, to speed next loop ? */ LWDEBUGF(1, "Edge %d is a known backward edge of the new ring", e->edge_id); backward_edges[backward_edges_count].edge_id = e->edge_id; backward_edges[backward_edges_count++].face_right = newface.face_id; found++; if ( found == 2 ) break; /* both edge sides are found on the ring */ } } if ( found ) continue; LWDEBUGF(1, "Edge %d is not a known edge of the new ring", e->edge_id); /* Check if the edge is now binding a different face */ if ( ! getPoint2d_p(e->geom->points, 0, &ep) ) { lwfree(signed_edge_ids); lwpoly_free(shell); lwfree(forward_edges); /* contents owned by edges */ lwfree(backward_edges); /* contents owned by edges */ _lwt_release_edges(edges, numfaceedges); lwerror("Edge %d is empty", e->edge_id); return -2; } /* IDEA: check that bounding box shortcut is taken, or use * shellbox to do it here */ contains = ptarray_contains_point(pa, &ep); LWDEBUGF(1, "Edge %d first point %s new ring", e->edge_id, (contains == LW_INSIDE ? "inside" : contains == LW_OUTSIDE ? "outside" : "on boundary of")); /* (2.2) skip edges (NOT, if newface_outside) contained in ring */ if ( newface_outside ) { if ( contains != LW_OUTSIDE ) { LWDEBUGF(1, "Edge %d not outside of the new ring, not updating it", e->edge_id); continue; } } else { if ( contains != LW_INSIDE ) { LWDEBUGF(1, "Edge %d not inside the new ring, not updating it", e->edge_id); continue; } } /* (2.3) push to forward_edges if left_face = oface */ if ( e->face_left == face ) { LWDEBUGF(1, "Edge %d has new face on the left side", e->edge_id); forward_edges[forward_edges_count].edge_id = e->edge_id; forward_edges[forward_edges_count++].face_left = newface.face_id; } /* (2.4) push to backward_edges if right_face = oface */ if ( e->face_right == face ) { LWDEBUGF(1, "Edge %d has new face on the right side", e->edge_id); backward_edges[backward_edges_count].edge_id = e->edge_id; backward_edges[backward_edges_count++].face_right = newface.face_id; } } /* 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( signed_edge_ids ); lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -2; } if ( ret != forward_edges_count ) { lwfree( signed_edge_ids ); lwerror("Unexpected error: %d edges updated when expecting %d", ret, forward_edges_count); return -2; } } /* 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( signed_edge_ids ); lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -2; } if ( ret != backward_edges_count ) { lwfree( signed_edge_ids ); lwerror("Unexpected error: %d edges updated when expecting %d", ret, backward_edges_count); return -2; } } lwfree(forward_edges); lwfree(backward_edges); } _lwt_release_edges(edges, numfaceedges); /* Update isolated nodes which are now in new face */ uint64_t numisonodes = 1; fields = LWT_COL_NODE_NODE_ID | LWT_COL_NODE_GEOM; LWT_ISO_NODE *nodes = lwt_be_getNodeByFace(topo, &face, &numisonodes, fields, newface.mbr); if (numisonodes == UINT64_MAX) { lwfree(signed_edge_ids); lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -2; } if ( numisonodes ) { LWT_ISO_NODE *updated_nodes = lwalloc(sizeof(LWT_ISO_NODE)*numisonodes); int nodes_to_update = 0; for (i=0; igeom->point, 0); int contains = ptarray_contains_point(pa, pt) == LW_INSIDE; LWDEBUGF(1, "Node %d is %scontained in new ring, newface is %s", n->node_id, contains ? "" : "not ", newface_outside ? "outside" : "inside" ); if ( newface_outside ) { if ( contains ) { LWDEBUGF(1, "Node %d contained in an hole of the new face", n->node_id); continue; } } else { if ( ! contains ) { LWDEBUGF(1, "Node %d not contained in the face shell", n->node_id); continue; } } updated_nodes[nodes_to_update].node_id = n->node_id; updated_nodes[nodes_to_update++].containing_face = newface.face_id; LWDEBUGF(1, "Node %d will be updated", n->node_id); } _lwt_release_nodes(nodes, numisonodes); if ( nodes_to_update ) { int ret = lwt_be_updateNodesById(topo, updated_nodes, nodes_to_update, LWT_COL_NODE_CONTAINING_FACE); if ( ret == -1 ) { lwfree( signed_edge_ids ); lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -2; } } lwfree(updated_nodes); } lwfree(signed_edge_ids); lwpoly_free(shell); 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, LWLINE *geom, int skipChecks, int modFace ) { LWT_ISO_EDGE newedge; LWGEOM *cleangeom; edgeend span; /* start point analisys */ edgeend epan; /* end point analisys */ POINT2D p1, pn, p2; POINTARRAY *pa; LWT_ELEMID node_ids[2]; const LWPOINT *start_node_geom = NULL; const LWPOINT *end_node_geom = NULL; uint64_t num_nodes; LWT_ISO_NODE *endpoints; uint64_t i; int prev_left; int prev_right; LWT_ISO_EDGE seledge; LWT_ISO_EDGE updedge; if ( ! skipChecks ) { /* curve must be simple */ if ( ! lwgeom_is_simple(lwline_as_lwgeom(geom)) ) { lwerror("SQL/MM Spatial exception - curve not simple"); return -1; } } newedge.start_node = start_node; newedge.end_node = end_node; newedge.geom = geom; newedge.face_left = -1; newedge.face_right = -1; /* TODO: should do the repeated points removal in 2D space */ cleangeom = lwgeom_remove_repeated_points( lwline_as_lwgeom(geom), 0 ); pa = lwgeom_as_lwline(cleangeom)->points; if ( pa->npoints < 2 ) { lwgeom_free(cleangeom); lwerror("Invalid edge (no two distinct vertices exist)"); return -1; } /* Initialize endpoint info (some of that ) */ span.cwFace = span.ccwFace = epan.cwFace = epan.ccwFace = -1; /* Compute azimuth of first edge end on start node */ getPoint2d_p(pa, 0, &p1); if ( ! _lwt_FirstDistinctVertex2D(pa, &p1, 0, 1, &pn) ) { lwgeom_free(cleangeom); lwerror("Invalid edge (no two distinct vertices exist)"); return -1; } if ( ! azimuth_pt_pt(&p1, &pn, &span.myaz) ) { lwgeom_free(cleangeom); lwerror("error computing azimuth of first edgeend [%.15g %.15g,%.15g %.15g]", p1.x, p1.y, pn.x, pn.y); return -1; } LWDEBUGF(1, "edge's start node is %g,%g", p1.x, p1.y); /* Compute azimuth of last edge end on end node */ getPoint2d_p(pa, pa->npoints-1, &p2); if ( ! _lwt_FirstDistinctVertex2D(pa, &p2, pa->npoints-1, -1, &pn) ) { lwgeom_free(cleangeom); /* This should never happen as we checked the edge while computing first edgend */ lwerror("Invalid clean edge (no two distinct vertices exist) - should not happen"); return -1; } lwgeom_free(cleangeom); if ( ! azimuth_pt_pt(&p2, &pn, &epan.myaz) ) { lwerror("error computing azimuth of last edgeend [%.15g %.15g,%.15g %.15g]", p2.x, p2.y, pn.x, pn.y); return -1; } LWDEBUGF(1, "edge's end node is %g,%g", p2.x, p2.y); /* * Check endpoints existence, match with Curve geometry * and get face information (if any) */ if ( start_node != end_node ) { num_nodes = 2; node_ids[0] = start_node; node_ids[1] = end_node; } else { num_nodes = 1; node_ids[0] = start_node; } endpoints = lwt_be_getNodeById( topo, node_ids, &num_nodes, LWT_COL_NODE_ALL ); if (num_nodes == UINT64_MAX) { lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } for ( i=0; icontaining_face != -1 ) { if ( newedge.face_left == -1 ) { newedge.face_left = newedge.face_right = node->containing_face; } else if ( newedge.face_left != node->containing_face ) { _lwt_release_nodes(endpoints, num_nodes); lwerror("SQL/MM Spatial exception - geometry crosses an edge" " (endnodes in faces %" LWTFMT_ELEMID " and %" LWTFMT_ELEMID ")", newedge.face_left, node->containing_face); } } LWDEBUGF(1, "Node %d, with geom %p (looking for %d and %d)", node->node_id, node->geom, start_node, end_node); if ( node->node_id == start_node ) { start_node_geom = node->geom; } if ( node->node_id == end_node ) { end_node_geom = node->geom; } } if ( ! skipChecks ) { if ( ! start_node_geom ) { if ( num_nodes ) _lwt_release_nodes(endpoints, num_nodes); lwerror("SQL/MM Spatial exception - non-existent node"); return -1; } else { pa = start_node_geom->point; getPoint2d_p(pa, 0, &pn); if ( ! P2D_SAME_STRICT(&pn, &p1) ) { if ( num_nodes ) _lwt_release_nodes(endpoints, num_nodes); lwerror("SQL/MM Spatial exception" " - start node not geometry start point." //" - start node not geometry start point (%g,%g != %g,%g).", pn.x, pn.y, p1.x, p1.y ); return -1; } } if ( ! end_node_geom ) { if ( num_nodes ) _lwt_release_nodes(endpoints, num_nodes); lwerror("SQL/MM Spatial exception - non-existent node"); return -1; } else { pa = end_node_geom->point; getPoint2d_p(pa, 0, &pn); if ( ! P2D_SAME_STRICT(&pn, &p2) ) { if ( num_nodes ) _lwt_release_nodes(endpoints, num_nodes); lwerror("SQL/MM Spatial exception" " - end node not geometry end point." //" - end node not geometry end point (%g,%g != %g,%g).", pn.x, pn.y, p2.x, p2.y ); return -1; } } if ( num_nodes ) _lwt_release_nodes(endpoints, num_nodes); if ( _lwt_CheckEdgeCrossing( topo, start_node, end_node, geom, 0 ) ) return -1; } /* ! skipChecks */ /* * All checks passed, time to prepare the new edge */ newedge.edge_id = lwt_be_getNextEdgeId( topo ); if ( newedge.edge_id == -1 ) { lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } /* Find adjacent edges to each endpoint */ int isclosed = start_node == end_node; int found; found = _lwt_FindAdjacentEdges( topo, start_node, &span, isclosed ? &epan : NULL, -1 ); if ( found ) { span.was_isolated = 0; newedge.next_right = span.nextCW ? span.nextCW : -newedge.edge_id; prev_left = span.nextCCW ? -span.nextCCW : newedge.edge_id; LWDEBUGF(1, "New edge %d is connected on start node, " "next_right is %d, prev_left is %d", newedge.edge_id, newedge.next_right, prev_left); if ( modFace != -1 ) { if ( newedge.face_right == -1 ) { newedge.face_right = span.cwFace; } if ( newedge.face_left == -1 ) { newedge.face_left = span.ccwFace; } } } else { span.was_isolated = 1; newedge.next_right = isclosed ? -newedge.edge_id : newedge.edge_id; prev_left = isclosed ? newedge.edge_id : -newedge.edge_id; LWDEBUGF(1, "New edge %d is isolated on start node, " "next_right is %d, prev_left is %d", newedge.edge_id, newedge.next_right, prev_left); } found = _lwt_FindAdjacentEdges( topo, end_node, &epan, isclosed ? &span : NULL, -1 ); if ( found ) { epan.was_isolated = 0; newedge.next_left = epan.nextCW ? epan.nextCW : newedge.edge_id; prev_right = epan.nextCCW ? -epan.nextCCW : -newedge.edge_id; LWDEBUGF(1, "New edge %d is connected on end node, " "next_left is %d, prev_right is %d", newedge.edge_id, newedge.next_left, prev_right); if ( modFace != -1 ) { if ( newedge.face_right == -1 ) { newedge.face_right = span.ccwFace; } else if ( newedge.face_right != epan.ccwFace ) { /* side-location conflict */ lwerror("Side-location conflict: " "new edge starts in face" " %" LWTFMT_ELEMID " and ends in face" " %" LWTFMT_ELEMID, newedge.face_right, epan.ccwFace ); return -1; } if ( newedge.face_left == -1 ) { newedge.face_left = span.cwFace; } else if ( newedge.face_left != epan.cwFace ) { /* side-location conflict */ lwerror("Side-location conflict: " "new edge starts in face" " %" LWTFMT_ELEMID " and ends in face" " %" LWTFMT_ELEMID, newedge.face_left, epan.cwFace ); return -1; } } } else { epan.was_isolated = 1; newedge.next_left = isclosed ? newedge.edge_id : -newedge.edge_id; prev_right = isclosed ? -newedge.edge_id : newedge.edge_id; LWDEBUGF(1, "New edge %d is isolated on end node, " "next_left is %d, prev_right is %d", newedge.edge_id, newedge.next_left, prev_right); } /* * If we don't have faces setup by now we must have encountered * a malformed topology (no containing_face on isolated nodes, no * left/right faces on adjacent edges or mismatching values) */ if ( modFace > -1 ) { if ( newedge.face_left != newedge.face_right ) { lwerror("Left(%" LWTFMT_ELEMID ")/right(%" LWTFMT_ELEMID ")" " faces mismatch: invalid topology ?", newedge.face_left, newedge.face_right); return -1; } else if ( newedge.face_left == -1 ) { lwerror("Could not derive edge face from linked primitives:" " invalid topology ?"); return -1; } } /* * Insert the new edge, and update all linking */ int ret = lwt_be_insertEdges(topo, &newedge, 1); if ( ret == -1 ) { lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } else if ( ret == 0 ) { lwerror("Insertion of split edge failed (no reason)"); return -1; } int updfields; /* Link prev_left to us * (if it's not us already) */ if ( llabs(prev_left) != newedge.edge_id ) { if ( prev_left > 0 ) { /* its next_left_edge is us */ updfields = LWT_COL_EDGE_NEXT_LEFT; updedge.next_left = newedge.edge_id; seledge.edge_id = prev_left; } else { /* its next_right_edge is us */ updfields = LWT_COL_EDGE_NEXT_RIGHT; updedge.next_right = newedge.edge_id; seledge.edge_id = -prev_left; } ret = lwt_be_updateEdges(topo, &seledge, LWT_COL_EDGE_EDGE_ID, &updedge, updfields, NULL, 0); if ( ret == -1 ) { lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } } /* Link prev_right to us * (if it's not us already) */ if ( llabs(prev_right) != newedge.edge_id ) { if ( prev_right > 0 ) { /* its next_left_edge is -us */ updfields = LWT_COL_EDGE_NEXT_LEFT; updedge.next_left = -newedge.edge_id; seledge.edge_id = prev_right; } else { /* its next_right_edge is -us */ updfields = LWT_COL_EDGE_NEXT_RIGHT; updedge.next_right = -newedge.edge_id; seledge.edge_id = -prev_right; } ret = lwt_be_updateEdges(topo, &seledge, LWT_COL_EDGE_EDGE_ID, &updedge, updfields, NULL, 0); if ( ret == -1 ) { lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } } /* NOT IN THE SPECS... * set containing_face = null for start_node and end_node * if they where isolated * */ LWT_ISO_NODE updnode, selnode; updnode.containing_face = -1; if ( span.was_isolated ) { selnode.node_id = start_node; ret = lwt_be_updateNodes(topo, &selnode, LWT_COL_NODE_NODE_ID, &updnode, LWT_COL_NODE_CONTAINING_FACE, NULL, 0); if ( ret == -1 ) { lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } } if ( epan.was_isolated ) { selnode.node_id = end_node; ret = lwt_be_updateNodes(topo, &selnode, LWT_COL_NODE_NODE_ID, &updnode, LWT_COL_NODE_CONTAINING_FACE, NULL, 0); if ( ret == -1 ) { lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } } /* 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"); return newedge.edge_id; /* no split */ } int newface1 = -1; /* IDEA: avoid building edge ring if input is closed, which means we * know in advance it splits a face */ if ( ! modFace ) { newface1 = _lwt_AddFaceSplit( topo, -newedge.edge_id, newedge.face_left, 0 ); if ( newface1 == 0 ) { LWDEBUG(1, "New edge does not split any face"); return newedge.edge_id; /* no split */ } } int newface = _lwt_AddFaceSplit( topo, newedge.edge_id, newedge.face_left, 0 ); if ( modFace ) { if ( newface == 0 ) { LWDEBUG(1, "New edge does not split any face"); return newedge.edge_id; /* no split */ } if ( newface < 0 ) { /* face on the left is the universe face */ /* must be forming a maximal ring in universal face */ newface = _lwt_AddFaceSplit( topo, -newedge.edge_id, newedge.face_left, 0 ); if ( newface < 0 ) return newedge.edge_id; /* no split */ } else { _lwt_AddFaceSplit( topo, -newedge.edge_id, newedge.face_left, 1 ); } } /* * Update topogeometries, if needed */ if ( newedge.face_left != 0 ) { ret = lwt_be_updateTopoGeomFaceSplit(topo, newedge.face_left, newface, newface1); if ( ret == 0 ) { lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } if ( ! modFace ) { /* drop old face from the face table */ ret = lwt_be_deleteFacesById(topo, &(newedge.face_left), 1); if ( ret == -1 ) { lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } } } } // end of face split checking return newedge.edge_id; } LWT_ELEMID lwt_AddEdgeModFace( LWT_TOPOLOGY* topo, LWT_ELEMID start_node, LWT_ELEMID end_node, LWLINE *geom, int skipChecks ) { return _lwt_AddEdge( topo, start_node, end_node, geom, skipChecks, 1 ); } LWT_ELEMID lwt_AddEdgeNewFaces( LWT_TOPOLOGY* topo, LWT_ELEMID start_node, LWT_ELEMID end_node, LWLINE *geom, int skipChecks ) { return _lwt_AddEdge( topo, start_node, end_node, geom, skipChecks, 0 ); } static LWGEOM * _lwt_FaceByEdges(LWT_TOPOLOGY *topo, LWT_ISO_EDGE *edges, int numfaceedges) { LWGEOM *outg; LWCOLLECTION *bounds; LWGEOM **geoms = lwalloc( sizeof(LWGEOM*) * numfaceedges ); int i, validedges = 0; for ( i=0; isrid, topo->hasZ, 0) ); } bounds = lwcollection_construct(MULTILINETYPE, topo->srid, NULL, /* gbox */ validedges, geoms); outg = lwgeom_buildarea( lwcollection_as_lwgeom(bounds) ); lwcollection_release(bounds); lwfree(geoms); #if 0 { size_t sz; char *wkt = lwgeom_to_wkt(outg, WKT_EXTENDED, 2, &sz); LWDEBUGF(1, "_lwt_FaceByEdges returning area: %s", wkt); lwfree(wkt); } #endif return outg; } LWGEOM* lwt_GetFaceGeometry(LWT_TOPOLOGY* topo, LWT_ELEMID faceid) { uint64_t numfaceedges; LWT_ISO_EDGE *edges; LWT_ISO_FACE *face; LWPOLY *out; LWGEOM *outg; uint64_t i, edgeid; int fields; if (faceid == 0) { lwerror("SQL/MM Spatial exception - universal face has no geometry"); return NULL; } /* Construct the face geometry */ numfaceedges = 1; fields = LWT_COL_EDGE_GEOM | LWT_COL_EDGE_EDGE_ID | LWT_COL_EDGE_FACE_LEFT | LWT_COL_EDGE_FACE_RIGHT ; edges = lwt_be_getEdgeByFace( topo, &faceid, &numfaceedges, fields, NULL ); if (numfaceedges == UINT64_MAX) { lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return NULL; } LWDEBUGF(1, "lwt_GetFaceGeometry: lwt_be_getEdgeByFace returned %d edges", numfaceedges); if ( numfaceedges == 0 ) { i = 1; face = lwt_be_getFaceById(topo, &faceid, &i, LWT_COL_FACE_FACE_ID); if (i == UINT64_MAX) { lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return NULL; } if ( i == 0 ) { lwerror("SQL/MM Spatial exception - non-existent face."); return NULL; } lwfree( face ); if ( i > 1 ) { lwerror("Corrupted topology: multiple face records have face_id=%" LWTFMT_ELEMID, faceid); return NULL; } /* Face has no boundary edges, we'll return EMPTY, see * https://trac.osgeo.org/postgis/ticket/3221 */ lwnotice("Corrupted topology: face %" LWTFMT_ELEMID " has no associated edges.", faceid); out = lwpoly_construct_empty(topo->srid, topo->hasZ, 0); return lwpoly_as_lwgeom(out); } edgeid = edges[0].edge_id; outg = _lwt_FaceByEdges( topo, edges, numfaceedges ); _lwt_release_edges(edges, numfaceedges); if ( ! outg ) { /* Face did have edges but no polygon could be constructed * with that material, sounds like a corrupted topology.. * * We'll return EMPTY, see * https://trac.osgeo.org/postgis/ticket/3221 */ lwnotice("Corrupted topology: face %" LWTFMT_ELEMID " could not be constructed only from edges " "knowing about it (like edge %" LWTFMT_ELEMID ").", faceid, edgeid); out = lwpoly_construct_empty(topo->srid, topo->hasZ, 0); return lwpoly_as_lwgeom(out); } return outg; } /* Find which edge from the "edges" set defines the next * portion of the given "ring". * * The edge might be either forward or backward. * * @param ring The ring to find definition of. * It is assumed it does not contain duplicated vertices. * @param from offset of the ring point to start looking from * @param edges array of edges to search into * @param numedges number of edges in the edges array * * @return index of the edge defining the next ring portion or * -1 if no edge was found to be part of the ring */ static int _lwt_FindNextRingEdge(const POINTARRAY *ring, int from, const LWT_ISO_EDGE *edges, int numedges) { int i; POINT2D p1; /* Get starting ring point */ getPoint2d_p(ring, from, &p1); LWDEBUGF(1, "Ring's 'from' point (%d) is %g,%g", from, p1.x, p1.y); /* find the edges defining the next portion of ring starting from * vertex "from" */ for ( i=0; igeom; POINTARRAY *epa = edge->points; POINT2D p2, pt; int match = 0; uint32_t j; /* Skip if the edge is a dangling one */ if ( isoe->face_left == isoe->face_right ) { LWDEBUGF(3, "_lwt_FindNextRingEdge: edge %" LWTFMT_ELEMID " has same face (%" LWTFMT_ELEMID ") on both sides, skipping", isoe->edge_id, isoe->face_left); continue; } if (epa->npoints < 2) { LWDEBUGF(3, "_lwt_FindNextRingEdge: edge %" LWTFMT_ELEMID " has only %"PRIu32" points", isoe->edge_id, epa->npoints); continue; } #if 0 size_t sz; LWDEBUGF(1, "Edge %" LWTFMT_ELEMID " is %s", isoe->edge_id, lwgeom_to_wkt(lwline_as_lwgeom(edge), WKT_EXTENDED, 2, &sz)); #endif /* ptarray_remove_repeated_points ? */ getPoint2d_p(epa, 0, &p2); LWDEBUGF(1, "Edge %" LWTFMT_ELEMID " 'first' point is %g,%g", isoe->edge_id, p2.x, p2.y); LWDEBUGF(1, "Rings's 'from' point is still %g,%g", p1.x, p1.y); if ( P2D_SAME_STRICT(&p1, &p2) ) { LWDEBUG(1, "P2D_SAME_STRICT(p1,p2) returned true"); LWDEBUGF(1, "First point of edge %" LWTFMT_ELEMID " matches ring vertex %d", isoe->edge_id, from); /* first point matches, let's check next non-equal one */ for ( j=1; jnpoints; ++j ) { getPoint2d_p(epa, j, &p2); LWDEBUGF(1, "Edge %" LWTFMT_ELEMID " 'next' point %d is %g,%g", isoe->edge_id, j, p2.x, p2.y); /* we won't check duplicated edge points */ if ( P2D_SAME_STRICT(&p1, &p2) ) continue; /* we assume there are no duplicated points in ring */ getPoint2d_p(ring, from+1, &pt); LWDEBUGF(1, "Ring's point %d is %g,%g", from+1, pt.x, pt.y); match = P2D_SAME_STRICT(&pt, &p2); break; /* we want to check a single non-equal next vertex */ } #if POSTGIS_DEBUG_LEVEL > 0 if ( match ) { LWDEBUGF(1, "Prev point of edge %" LWTFMT_ELEMID " matches ring vertex %d", isoe->edge_id, from+1); } else { LWDEBUGF(1, "Prev point of edge %" LWTFMT_ELEMID " does not match ring vertex %d", isoe->edge_id, from+1); } #endif } if ( ! match ) { LWDEBUGF(1, "Edge %" LWTFMT_ELEMID " did not match as forward", isoe->edge_id); getPoint2d_p(epa, epa->npoints-1, &p2); LWDEBUGF(1, "Edge %" LWTFMT_ELEMID " 'last' point is %g,%g", isoe->edge_id, p2.x, p2.y); if ( P2D_SAME_STRICT(&p1, &p2) ) { LWDEBUGF(1, "Last point of edge %" LWTFMT_ELEMID " matches ring vertex %d", isoe->edge_id, from); /* last point matches, let's check next non-equal one */ for ( j=2; j<=epa->npoints; j++ ) { getPoint2d_p(epa, epa->npoints - j, &p2); LWDEBUGF(1, "Edge %" LWTFMT_ELEMID " 'prev' point %d is %g,%g", isoe->edge_id, epa->npoints - j, p2.x, p2.y); /* we won't check duplicated edge points */ if ( P2D_SAME_STRICT(&p1, &p2) ) continue; /* we assume there are no duplicated points in ring */ getPoint2d_p(ring, from+1, &pt); LWDEBUGF(1, "Ring's point %d is %g,%g", from+1, pt.x, pt.y); match = P2D_SAME_STRICT(&pt, &p2); break; /* we want to check a single non-equal next vertex */ } } #if POSTGIS_DEBUG_LEVEL > 0 if ( match ) { LWDEBUGF(1, "Prev point of edge %" LWTFMT_ELEMID " matches ring vertex %d", isoe->edge_id, from+1); } else { LWDEBUGF(1, "Prev point of edge %" LWTFMT_ELEMID " does not match ring vertex %d", isoe->edge_id, from+1); } #endif } if ( match ) return i; } return -1; } /* Reverse values in array between "from" (inclusive) * and "to" (exclusive) indexes */ static void _lwt_ReverseElemidArray(LWT_ELEMID *ary, int from, int to) { LWT_ELEMID t; while (from < to) { t = ary[from]; ary[from++] = ary[to]; ary[to--] = t; } } /* Rotate values in array between "from" (inclusive) * and "to" (exclusive) indexes, so that "rotidx" is * the new value at "from" */ static void _lwt_RotateElemidArray(LWT_ELEMID *ary, int from, int to, int rotidx) { _lwt_ReverseElemidArray(ary, from, rotidx-1); _lwt_ReverseElemidArray(ary, rotidx, to-1); _lwt_ReverseElemidArray(ary, from, to-1); } int lwt_GetFaceEdges(LWT_TOPOLOGY* topo, LWT_ELEMID face_id, LWT_ELEMID **out ) { LWGEOM *face; LWPOLY *facepoly; LWT_ISO_EDGE *edges; uint64_t numfaceedges; int fields; uint32_t i; int nseid = 0; /* number of signed edge ids */ int prevseid; LWT_ELEMID *seid; /* signed edge ids */ /* Get list of face edges */ numfaceedges = 1; fields = LWT_COL_EDGE_EDGE_ID | LWT_COL_EDGE_GEOM | LWT_COL_EDGE_FACE_LEFT | LWT_COL_EDGE_FACE_RIGHT ; edges = lwt_be_getEdgeByFace( topo, &face_id, &numfaceedges, fields, NULL ); if (numfaceedges == UINT64_MAX) { lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } if ( ! numfaceedges ) return 0; /* no edges in output */ LWDEBUGF(1, "lwt_GetFaceEdges: lwt_be_getEdgeByFace returned %d edges", numfaceedges); /* order edges by occurrence in face */ face = _lwt_FaceByEdges(topo, edges, numfaceedges); if ( ! face ) { /* _lwt_FaceByEdges should have already invoked lwerror in this case */ _lwt_release_edges(edges, numfaceedges); return -1; } if ( lwgeom_is_empty(face) ) { /* no edges in output */ _lwt_release_edges(edges, numfaceedges); lwgeom_free(face); return 0; } /* force_lhr, if the face is not the universe */ /* _lwt_FaceByEdges seems to guaranteed RHR */ /* lwgeom_force_clockwise(face); */ if ( face_id ) lwgeom_reverse_in_place(face); #if 0 { size_t sz; char *wkt = lwgeom_to_wkt(face, WKT_EXTENDED, 6, &sz); LWDEBUGF(1, "Geometry of face %" LWTFMT_ELEMID " is: %s", face_id, wkt); lwfree(wkt); } #endif facepoly = lwgeom_as_lwpoly(face); if ( ! facepoly ) { _lwt_release_edges(edges, numfaceedges); lwgeom_free(face); lwerror("Geometry of face %" LWTFMT_ELEMID " is not a polygon", face_id); return -1; } nseid = prevseid = 0; seid = lwalloc( sizeof(LWT_ELEMID) * numfaceedges ); /* for each ring of the face polygon... */ for ( i=0; inrings; ++i ) { const POINTARRAY *ring = facepoly->rings[i]; int32_t j = 0; LWT_ISO_EDGE *nextedge; LWLINE *nextline; LWDEBUGF(1, "Ring %d has %d points", i, ring->npoints); while ( j < (int32_t) ring->npoints-1 ) { LWDEBUGF(1, "Looking for edge covering ring %d from vertex %d", i, j); int edgeno = _lwt_FindNextRingEdge(ring, j, edges, numfaceedges); if ( edgeno == -1 ) { /* should never happen */ _lwt_release_edges(edges, numfaceedges); lwgeom_free(face); lwfree(seid); lwerror("No edge (among %d) found to be defining geometry of face %" LWTFMT_ELEMID, numfaceedges, face_id); return -1; } nextedge = &(edges[edgeno]); nextline = nextedge->geom; LWDEBUGF(1, "Edge %" LWTFMT_ELEMID " covers ring %d from vertex %d to %d", nextedge->edge_id, i, j, j + nextline->points->npoints - 1); #if 0 { size_t sz; char *wkt = lwgeom_to_wkt(lwline_as_lwgeom(nextline), WKT_EXTENDED, 6, &sz); LWDEBUGF(1, "Edge %" LWTFMT_ELEMID " is %s", nextedge->edge_id, wkt); lwfree(wkt); } #endif j += nextline->points->npoints - 1; /* Add next edge to the output array */ seid[nseid++] = nextedge->face_left == face_id ? nextedge->edge_id : -nextedge->edge_id; /* avoid checking again on next time turn */ nextedge->face_left = nextedge->face_right = -1; } /* now "scroll" the list of edges so that the one * with smaller absolute edge_id is first */ /* Range is: [prevseid, nseid) -- [inclusive, exclusive) */ if ( (nseid - prevseid) > 1 ) {{ LWT_ELEMID minid = 0; int minidx = 0; LWDEBUGF(1, "Looking for smallest id among the %d edges " "composing ring %d", (nseid-prevseid), i); for ( j=prevseid; jbe_iface)); return -1; } else if ( i == 0 ) { lwerror("SQL/MM Spatial exception - non-existent edge %" LWTFMT_ELEMID, edge_id); return -1; } else { lwerror("Backend coding error: getEdgeById callback returned NULL " "but numelements output parameter has value %d " "(expected 0 or 1)", i); return -1; } } LWDEBUGF(1, "lwt_ChangeEdgeGeom: " "old edge has %d points, new edge has %d points", oldedge->geom->points->npoints, geom->points->npoints); /* * e) Check StartPoint consistency */ getPoint2d_p(oldedge->geom->points, 0, &p1); getPoint2d_p(geom->points, 0, &pt); if ( ! P2D_SAME_STRICT(&p1, &pt) ) { _lwt_release_edges(oldedge, 1); lwerror("SQL/MM Spatial exception - " "start node not geometry start point."); return -1; } /* * f) Check EndPoint consistency */ if ( oldedge->geom->points->npoints < 2 ) { _lwt_release_edges(oldedge, 1); lwerror("Corrupted topology: edge %" LWTFMT_ELEMID " has less than 2 vertices", oldedge->edge_id); return -1; } getPoint2d_p(oldedge->geom->points, oldedge->geom->points->npoints-1, &p2); if ( geom->points->npoints < 2 ) { _lwt_release_edges(oldedge, 1); lwerror("Invalid edge: less than 2 vertices"); return -1; } getPoint2d_p(geom->points, geom->points->npoints-1, &pt); if ( ! P2D_SAME_STRICT(&pt, &p2) ) { _lwt_release_edges(oldedge, 1); lwerror("SQL/MM Spatial exception - " "end node not geometry end point."); return -1; } /* Not in the specs: * if the edge is closed, check we didn't change winding ! * (should be part of isomorphism checking) */ if ( oldedge->start_node == oldedge->end_node ) { isclosed = 1; #if 1 /* TODO: this is actually bogus as a test */ /* check for valid edge (distinct vertices must exist) */ if ( ! _lwt_GetInteriorEdgePoint(geom, &pt) ) { _lwt_release_edges(oldedge, 1); lwerror("Invalid edge (no two distinct vertices exist)"); return -1; } #endif if ( ptarray_isccw(oldedge->geom->points) != ptarray_isccw(geom->points) ) { _lwt_release_edges(oldedge, 1); lwerror("Edge twist at node POINT(%g %g)", p1.x, p1.y); return -1; } } if ( _lwt_CheckEdgeCrossing(topo, oldedge->start_node, oldedge->end_node, geom, edge_id ) ) { /* would have called lwerror already, leaking :( */ _lwt_release_edges(oldedge, 1); return -1; } LWDEBUG(1, "lwt_ChangeEdgeGeom: " "edge crossing check passed "); /* * Not in the specs: * Check topological isomorphism */ /* Check that the "motion range" doesn't include any node */ // 1. compute combined bbox of old and new edge GBOX mbox; /* motion box */ lwgeom_add_bbox((LWGEOM*)oldedge->geom); /* just in case */ lwgeom_add_bbox((LWGEOM*)geom); /* just in case */ gbox_union(oldedge->geom->bbox, geom->bbox, &mbox); // 2. fetch all nodes in the combined box LWT_ISO_NODE *nodes; uint64_t numnodes; nodes = lwt_be_getNodeWithinBox2D(topo, &mbox, &numnodes, LWT_COL_NODE_ALL, 0); LWDEBUGF(1, "lwt_be_getNodeWithinBox2D returned %d nodes", numnodes); if (numnodes == UINT64_MAX) { _lwt_release_edges(oldedge, 1); lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } // 3. if any node beside endnodes are found: if ( numnodes > ( 1 + isclosed ? 0 : 1 ) ) {{ // 3.2. bail out if any node is in one and not the other for (i=0; inode_id == oldedge->start_node ) continue; if ( n->node_id == oldedge->end_node ) continue; const POINT2D *pt = getPoint2d_cp(n->geom->point, 0); int ocont = ptarray_contains_point_partial(oldedge->geom->points, pt, isclosed, NULL) == LW_INSIDE; int ncont = ptarray_contains_point_partial(geom->points, pt, isclosed, NULL) == LW_INSIDE; if (ocont != ncont) { size_t sz; char *wkt = lwgeom_to_wkt(lwpoint_as_lwgeom(n->geom), WKT_ISO, 15, &sz); _lwt_release_nodes(nodes, numnodes); lwerror("Edge motion collision at %s", wkt); lwfree(wkt); /* would not necessarely reach this point */ return -1; } } }} if ( numnodes ) _lwt_release_nodes(nodes, numnodes); LWDEBUG(1, "nodes containment check passed"); /* * Check edge adjacency before * TODO: can be optimized to gather azimuths of all edge ends once */ edgeend span_pre, epan_pre; /* initialize span_pre.myaz and epan_pre.myaz with existing edge */ int res = _lwt_InitEdgeEndByLine(&span_pre, &epan_pre, oldedge->geom, &p1, &p2); if (res) return -1; /* lwerror should have been raised */ _lwt_FindAdjacentEdges( topo, oldedge->start_node, &span_pre, isclosed ? &epan_pre : NULL, edge_id ); _lwt_FindAdjacentEdges( topo, oldedge->end_node, &epan_pre, isclosed ? &span_pre : NULL, edge_id ); LWDEBUGF(1, "edges adjacent to old edge are %" LWTFMT_ELEMID " and %" LWTFMT_ELEMID " (first point), %" LWTFMT_ELEMID " and %" LWTFMT_ELEMID " (last point)", span_pre.nextCW, span_pre.nextCCW, epan_pre.nextCW, epan_pre.nextCCW); /* update edge geometry */ newedge.edge_id = edge_id; newedge.geom = geom; res = lwt_be_updateEdgesById(topo, &newedge, 1, LWT_COL_EDGE_GEOM); if (res == -1) { _lwt_release_edges(oldedge, 1); lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } if (!res) { _lwt_release_edges(oldedge, 1); lwerror("Unexpected error: %d edges updated when expecting 1", i); return -1; } /* * Check edge adjacency after */ edgeend span_post, epan_post; /* initialize epan_post.myaz and epan_post.myaz */ res = _lwt_InitEdgeEndByLine(&span_post, &epan_post, geom, &p1, &p2); if (res) return -1; /* lwerror should have been raised */ _lwt_FindAdjacentEdges( topo, oldedge->start_node, &span_post, isclosed ? &epan_post : NULL, edge_id ); _lwt_FindAdjacentEdges( topo, oldedge->end_node, &epan_post, isclosed ? &span_post : NULL, edge_id ); LWDEBUGF(1, "edges adjacent to new edge are %" LWTFMT_ELEMID " and %" LWTFMT_ELEMID " (first point), %" LWTFMT_ELEMID " and %" LWTFMT_ELEMID " (last point)", span_pre.nextCW, span_pre.nextCCW, epan_pre.nextCW, epan_pre.nextCCW); /* Bail out if next CW or CCW edge on start node changed */ if ( span_pre.nextCW != span_post.nextCW || span_pre.nextCCW != span_post.nextCCW ) {{ LWT_ELEMID nid = oldedge->start_node; _lwt_release_edges(oldedge, 1); lwerror("Edge changed disposition around start node %" LWTFMT_ELEMID, nid); return -1; }} /* Bail out if next CW or CCW edge on end node changed */ if ( epan_pre.nextCW != epan_post.nextCW || epan_pre.nextCCW != epan_post.nextCCW ) {{ LWT_ELEMID nid = oldedge->end_node; _lwt_release_edges(oldedge, 1); lwerror("Edge changed disposition around end node %" LWTFMT_ELEMID, nid); return -1; }} /* -- Update faces MBR of left and right faces -- TODO: think about ways to optimize this part, like see if -- the old edge geometry participated in the definition -- of the current MBR (for shrinking) or the new edge MBR -- would be larger than the old face MBR... -- */ LWGEOM *oldgeom = lwline_as_lwgeom(oldedge->geom); LWGEOM *newgeom = lwline_as_lwgeom(geom); lwgeom_refresh_bbox(oldgeom); /* Ensure we use a fit mbr, see #5709 -- TODO: fix this at lower level */ lwgeom_refresh_bbox(newgeom); /* Ensure we use a fit mbr, see #5709 -- TODO: fix this at lower level */ const GBOX* oldbox = lwgeom_get_bbox(oldgeom); const GBOX* newbox = lwgeom_get_bbox(newgeom); if ( ! gbox_same(oldbox, newbox) ) { GBOX* updatedBox; uint64_t facestoupdate = 0; LWT_ISO_FACE faces[2]; if ( oldedge->face_left > 0 ) { updatedBox = lwt_be_computeFaceMBR(topo, oldedge->face_left); if ( ! updatedBox ) { lwerror("Corrupted topology: face %d, left of edge %d, has no bbox", oldedge->face_left, edge_id); return -1; } faces[facestoupdate].face_id = oldedge->face_left; /* ownership transferred to faces[] */ faces[facestoupdate++].mbr = updatedBox; } if ( oldedge->face_right > 0 /* no need to update twice the same face.. */ && oldedge->face_right != oldedge->face_left ) { updatedBox = lwt_be_computeFaceMBR(topo, oldedge->face_right); if ( ! updatedBox ) { lwerror("Corrupted topology: face %d, right of edge %d, has no bbox", oldedge->face_right, edge_id); return -1; } faces[facestoupdate].face_id = oldedge->face_right; /* ownership transferred to faces[] */ faces[facestoupdate++].mbr = updatedBox; } LWDEBUGF(1, "%d faces to update", facestoupdate); if ( facestoupdate ) { uint64_t updatedFaces = lwt_be_updateFacesById(topo, &(faces[0]), facestoupdate); if (updatedFaces != facestoupdate) { while ( facestoupdate-- ) lwfree(faces[facestoupdate].mbr); _lwt_release_edges(oldedge, 1); if (updatedFaces == UINT64_MAX) lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); else lwerror("Unexpected error: %d faces updated when expecting 1", updatedFaces); return -1; } } while ( facestoupdate-- ) lwfree(faces[facestoupdate].mbr); } else { LWDEBUG(1, "BBOX of changed edge did not change"); } LWDEBUG(1, "all done, cleaning up edges"); _lwt_release_edges(oldedge, 1); return 0; /* success */ } /* Only return CONTAINING_FACE in the node object */ static LWT_ISO_NODE * _lwt_GetIsoNode(LWT_TOPOLOGY* topo, LWT_ELEMID nid) { LWT_ISO_NODE *node; uint64_t n = 1; node = lwt_be_getNodeById( topo, &nid, &n, LWT_COL_NODE_CONTAINING_FACE ); if (n == UINT64_MAX) { lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return 0; } if ( n < 1 ) { lwerror("SQL/MM Spatial exception - non-existent node"); return 0; } if ( node->containing_face == -1 ) { lwfree(node); lwerror("SQL/MM Spatial exception - not isolated node"); return 0; } return node; } int lwt_MoveIsoNode(LWT_TOPOLOGY* topo, LWT_ELEMID nid, LWPOINT *pt) { LWT_ISO_NODE *node; int ret; int newPointFace; node = _lwt_GetIsoNode( topo, nid ); if ( ! node ) return -1; if ( lwt_be_ExistsCoincidentNode(topo, pt) ) { lwfree(node); lwerror("SQL/MM Spatial exception - coincident node"); return -1; } if ( lwt_be_ExistsEdgeIntersectingPoint(topo, pt) ) { lwfree(node); lwerror("SQL/MM Spatial exception - edge crosses node."); return -1; } /* Check that the new point is in the same containing face ! * See https://trac.osgeo.org/postgis/ticket/3232 */ newPointFace = lwt_GetFaceContainingPoint(topo, pt); if ( newPointFace == -1 ) { lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } if ( node->containing_face != newPointFace ) { lwfree(node); lwerror("Cannot move isolated node across faces"); return -1; } node->node_id = nid; node->geom = pt; ret = lwt_be_updateNodesById(topo, node, 1, LWT_COL_NODE_GEOM); if ( ret == -1 ) { lwfree(node); lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } lwfree(node); return 0; } int lwt_RemoveIsoNode(LWT_TOPOLOGY* topo, LWT_ELEMID nid) { LWT_ISO_NODE *node; int n = 1; node = _lwt_GetIsoNode( topo, nid ); if ( ! node ) return -1; n = lwt_be_deleteNodesById( topo, &nid, n ); if ( n == -1 ) { lwfree(node); lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } if ( n != 1 ) { lwfree(node); lwerror("Unexpected error: %d nodes deleted when expecting 1", n); return -1; } if ( ! lwt_be_checkTopoGeomRemIsoNode(topo, nid) ) { lwfree(node); lwerror("%s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } lwfree(node); return 0; /* success */ } int lwt_RemIsoEdge(LWT_TOPOLOGY* topo, LWT_ELEMID id) { LWT_ISO_EDGE deledge; LWT_ISO_EDGE *edge; LWT_ELEMID nid[2]; LWT_ISO_NODE upd_node[2]; LWT_ELEMID containing_face; uint64_t n = 1; uint64_t i; edge = lwt_be_getEdgeById( topo, &id, &n, LWT_COL_EDGE_START_NODE| LWT_COL_EDGE_END_NODE | LWT_COL_EDGE_FACE_LEFT | LWT_COL_EDGE_FACE_RIGHT ); if ( ! edge ) { lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } if ( ! n ) { lwerror("SQL/MM Spatial exception - non-existent edge"); return -1; } if ( n > 1 ) { lwfree(edge); lwerror("Corrupted topology: more than a single edge have id %" LWTFMT_ELEMID, id); return -1; } if ( edge[0].face_left != edge[0].face_right ) { lwfree(edge); lwerror("SQL/MM Spatial exception - not isolated edge"); return -1; } containing_face = edge[0].face_left; nid[0] = edge[0].start_node; nid[1] = edge[0].end_node; lwfree(edge); n = 2; edge = lwt_be_getEdgeByNode( topo, nid, &n, LWT_COL_EDGE_EDGE_ID ); if ((n == UINT64_MAX) || (edge == NULL)) { lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } for (i = 0; i < n; ++i) { if (edge[i].edge_id != id) { lwfree(edge); lwerror("SQL/MM Spatial exception - not isolated edge"); return -1; } } lwfree(edge); deledge.edge_id = id; n = lwt_be_deleteEdges( topo, &deledge, LWT_COL_EDGE_EDGE_ID ); if (n == UINT64_MAX) { lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } if ( n != 1 ) { lwerror("Unexpected error: %d edges deleted when expecting 1", n); return -1; } upd_node[0].node_id = nid[0]; upd_node[0].containing_face = containing_face; n = 1; if ( nid[1] != nid[0] ) { upd_node[1].node_id = nid[1]; upd_node[1].containing_face = containing_face; ++n; } n = lwt_be_updateNodesById(topo, upd_node, n, LWT_COL_NODE_CONTAINING_FACE); if (n == UINT64_MAX) { lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } /* Check that the edge can be safely removed * See https://trac.osgeo.org/postgis/ticket/3248 */ if ( ! lwt_be_checkTopoGeomRemIsoEdge(topo, id) ) { lwerror("%s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } return 0; /* success */ } /* Used by _lwt_RemEdge to update edge face ref on healing * * @param of old face id (never 0 as you cannot remove face 0) * @param nf new face id * @return 0 on success, -1 on backend error */ static int _lwt_UpdateEdgeFaceRef( LWT_TOPOLOGY *topo, LWT_ELEMID of, LWT_ELEMID nf) { LWT_ISO_EDGE sel_edge, upd_edge; int ret; assert( of != 0 ); /* Update face_left for all edges still referencing old face */ sel_edge.face_left = of; upd_edge.face_left = nf; ret = lwt_be_updateEdges(topo, &sel_edge, LWT_COL_EDGE_FACE_LEFT, &upd_edge, LWT_COL_EDGE_FACE_LEFT, NULL, 0); if ( ret == -1 ) return -1; /* Update face_right for all edges still referencing old face */ sel_edge.face_right = of; upd_edge.face_right = nf; ret = lwt_be_updateEdges(topo, &sel_edge, LWT_COL_EDGE_FACE_RIGHT, &upd_edge, LWT_COL_EDGE_FACE_RIGHT, NULL, 0); if ( ret == -1 ) return -1; return 0; } /* Used by _lwt_RemEdge to update node face ref on healing * * @param of old face id (never 0 as you cannot remove face 0) * @param nf new face id * @return 0 on success, -1 on backend error */ static int _lwt_UpdateNodeFaceRef( LWT_TOPOLOGY *topo, LWT_ELEMID of, LWT_ELEMID nf) { LWT_ISO_NODE sel, upd; int ret; assert( of != 0 ); /* Update face_left for all edges still referencing old face */ sel.containing_face = of; upd.containing_face = nf; ret = lwt_be_updateNodes(topo, &sel, LWT_COL_NODE_CONTAINING_FACE, &upd, LWT_COL_NODE_CONTAINING_FACE, NULL, 0); if ( ret == -1 ) return -1; return 0; } /* Used by lwt_RemEdgeModFace and lwt_RemEdgeNewFaces * * Returns -1 on error, identifier of the face that takes up the space * previously occupied by the removed edge if modFace is 1, identifier of * the created face (0 if none) if modFace is 0. */ static LWT_ELEMID _lwt_RemEdge( LWT_TOPOLOGY* topo, LWT_ELEMID edge_id, int modFace ) { uint64_t i, nedges, nfaces, fields; LWT_ISO_EDGE *edge = NULL; LWT_ISO_EDGE *upd_edge = NULL; LWT_ISO_EDGE upd_edge_left[2]; int nedge_left = 0; LWT_ISO_EDGE upd_edge_right[2]; int nedge_right = 0; LWT_ISO_NODE upd_node[2]; int nnode = 0; LWT_ISO_FACE *faces = NULL; LWT_ISO_FACE newface; LWT_ELEMID node_ids[2]; LWT_ELEMID face_ids[2]; int fnode_edges = 0; /* number of edges on the first node (excluded * the one being removed ) */ int lnode_edges = 0; /* number of edges on the last node (excluded * the one being removed ) */ newface.face_id = 0; i = 1; edge = lwt_be_getEdgeById(topo, &edge_id, &i, LWT_COL_EDGE_ALL); if (!edge) { LWDEBUGF(1, "lwt_be_getEdgeById returned NULL and set i=%d", i); if (i == UINT64_MAX) { lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } else if (i == 0) { lwerror("SQL/MM Spatial exception - non-existent edge %" LWTFMT_ELEMID, edge_id); return -1; } else { lwerror( "Backend coding error: getEdgeById callback returned NULL " "but numelements output parameter has value %d " "(expected 0 or 1)", i); return -1; } } if ( ! lwt_be_checkTopoGeomRemEdge(topo, edge_id, edge->face_left, edge->face_right) ) { lwerror("%s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } LWDEBUG(1, "Updating next_{right,left}_face of ring edges..."); /* Update edge linking */ nedges = 0; node_ids[nedges++] = edge->start_node; if ( edge->end_node != edge->start_node ) { node_ids[nedges++] = edge->end_node; } fields = LWT_COL_EDGE_EDGE_ID | LWT_COL_EDGE_START_NODE | LWT_COL_EDGE_END_NODE | LWT_COL_EDGE_NEXT_LEFT | LWT_COL_EDGE_NEXT_RIGHT; upd_edge = lwt_be_getEdgeByNode( topo, &(node_ids[0]), &nedges, fields ); if (nedges == UINT64_MAX) { lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } nedge_left = nedge_right = 0; for ( i=0; iedge_id == edge_id ) continue; if ( e->start_node == edge->start_node || e->end_node == edge->start_node ) { ++fnode_edges; } if ( e->start_node == edge->end_node || e->end_node == edge->end_node ) { ++lnode_edges; } if ( e->next_left == -edge_id ) { upd_edge_left[nedge_left].edge_id = e->edge_id; upd_edge_left[nedge_left++].next_left = edge->next_left != edge_id ? edge->next_left : edge->next_right; } else if ( e->next_left == edge_id ) { upd_edge_left[nedge_left].edge_id = e->edge_id; upd_edge_left[nedge_left++].next_left = edge->next_right != -edge_id ? edge->next_right : edge->next_left; } if ( e->next_right == -edge_id ) { upd_edge_right[nedge_right].edge_id = e->edge_id; upd_edge_right[nedge_right++].next_right = edge->next_left != edge_id ? edge->next_left : edge->next_right; } else if ( e->next_right == edge_id ) { upd_edge_right[nedge_right].edge_id = e->edge_id; upd_edge_right[nedge_right++].next_right = edge->next_right != -edge_id ? edge->next_right : edge->next_left; } } if ( nedge_left ) { LWDEBUGF(1, "updating %d 'next_left' edges", nedge_left); /* update edges in upd_edge_left set next_left */ int result = lwt_be_updateEdgesById(topo, &(upd_edge_left[0]), nedge_left, LWT_COL_EDGE_NEXT_LEFT); if (result == -1) { _lwt_release_edges(edge, 1); lwfree(upd_edge); lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } } if ( nedge_right ) { LWDEBUGF(1, "updating %d 'next_right' edges", nedge_right); /* update edges in upd_edge_right set next_right */ int result = lwt_be_updateEdgesById(topo, &(upd_edge_right[0]), nedge_right, LWT_COL_EDGE_NEXT_RIGHT); if (result == -1) { _lwt_release_edges(edge, 1); lwfree(upd_edge); lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } } LWDEBUGF(1, "releasing %d updateable edges in %p", nedges, upd_edge); lwfree(upd_edge); /* Id of face that will take up all the space previously * taken by left and right faces of the edge */ LWT_ELEMID floodface; /* Find floodface, and update its mbr if != 0 */ if ( edge->face_left == edge->face_right ) { floodface = edge->face_right; } else { /* Two faces healed */ if ( edge->face_left == 0 || edge->face_right == 0 ) { floodface = 0; LWDEBUG(1, "floodface is universe"); } else { /* we choose right face as the face that will remain * to be symmetric with ST_AddEdgeModFace */ floodface = edge->face_right; LWDEBUGF(1, "floodface is %" LWTFMT_ELEMID, floodface); /* update mbr of floodface as union of mbr of both faces */ face_ids[0] = edge->face_left; face_ids[1] = edge->face_right; nfaces = 2; fields = LWT_COL_FACE_ALL; faces = lwt_be_getFaceById(topo, face_ids, &nfaces, fields); if (nfaces == UINT64_MAX) { lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } GBOX *box1=NULL; GBOX *box2=NULL; for ( i=0; iface_left ) { if ( ! box1 ) box1 = faces[i].mbr; else { i = edge->face_left; _lwt_release_edges(edge, 1); _lwt_release_faces(faces, nfaces); lwerror("corrupted topology: more than 1 face have face_id=%" LWTFMT_ELEMID, i); return -1; } } else if ( faces[i].face_id == edge->face_right ) { if ( ! box2 ) box2 = faces[i].mbr; else { i = edge->face_right; _lwt_release_edges(edge, 1); _lwt_release_faces(faces, nfaces); lwerror("corrupted topology: more than 1 face have face_id=%" LWTFMT_ELEMID, i); return -1; } } else { i = faces[i].face_id; _lwt_release_edges(edge, 1); _lwt_release_faces(faces, nfaces); lwerror("Backend coding error: getFaceById returned face " "with non-requested id %" LWTFMT_ELEMID, i); return -1; } } if ( ! box1 ) { i = edge->face_left; _lwt_release_edges(edge, 1); if ( nfaces ) _lwt_release_faces(faces, nfaces); lwerror("corrupted topology: no face have face_id=%" LWTFMT_ELEMID " (left face for edge %" LWTFMT_ELEMID ")", i, edge_id); return -1; } if ( ! box2 ) { i = edge->face_right; _lwt_release_edges(edge, 1); if ( nfaces ) _lwt_release_faces(faces, nfaces); lwerror("corrupted topology: no face have face_id=%" LWTFMT_ELEMID " (right face for edge %" LWTFMT_ELEMID ")", i, edge_id); return -1; } gbox_merge(box2, box1); /* box1 is now the union of the two */ newface.mbr = box1; if ( modFace ) { newface.face_id = floodface; int result = lwt_be_updateFacesById(topo, &newface, 1); _lwt_release_faces(faces, 2); if (result == -1) { _lwt_release_edges(edge, 1); lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } if (result != 1) { _lwt_release_edges(edge, 1); lwerror("Unexpected error: %d faces updated when expecting 1", i); return -1; } } else { /* New face replaces the old two faces */ newface.face_id = -1; int result = lwt_be_insertFaces(topo, &newface, 1); _lwt_release_faces(faces, 2); if (result == -1) { _lwt_release_edges(edge, 1); lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } if (result != 1) { _lwt_release_edges(edge, 1); lwerror("Unexpected error: %d faces inserted when expecting 1", result); return -1; } floodface = newface.face_id; } } /* Update face references for edges and nodes still referencing * the removed face(s) */ if ( edge->face_left != floodface ) { if ( -1 == _lwt_UpdateEdgeFaceRef(topo, edge->face_left, floodface) ) { _lwt_release_edges(edge, 1); lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } if ( -1 == _lwt_UpdateNodeFaceRef(topo, edge->face_left, floodface) ) { _lwt_release_edges(edge, 1); lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } } if ( edge->face_right != floodface ) { if ( -1 == _lwt_UpdateEdgeFaceRef(topo, edge->face_right, floodface) ) { _lwt_release_edges(edge, 1); lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } if ( -1 == _lwt_UpdateNodeFaceRef(topo, edge->face_right, floodface) ) { _lwt_release_edges(edge, 1); lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } } /* Update topogeoms on heal */ if ( ! lwt_be_updateTopoGeomFaceHeal(topo, edge->face_right, edge->face_left, floodface) ) { _lwt_release_edges(edge, 1); lwerror("%s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } } /* two faces healed */ /* Delete the edge */ int result = lwt_be_deleteEdges(topo, edge, LWT_COL_EDGE_EDGE_ID); if (result == -1) { _lwt_release_edges(edge, 1); lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } /* If any of the edge nodes remained isolated, set * containing_face = floodface */ if ( ! fnode_edges ) { upd_node[nnode].node_id = edge->start_node; upd_node[nnode].containing_face = floodface; ++nnode; } if ( edge->end_node != edge->start_node && ! lnode_edges ) { upd_node[nnode].node_id = edge->end_node; upd_node[nnode].containing_face = floodface; ++nnode; } if ( nnode ) { int result = lwt_be_updateNodesById(topo, upd_node, nnode, LWT_COL_NODE_CONTAINING_FACE); if (result == -1) { _lwt_release_edges(edge, 1); lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } } if ( edge->face_left != edge->face_right ) /* or there'd be no face to remove */ { LWT_ELEMID ids[2]; int nids = 0; if ( edge->face_right != floodface ) ids[nids++] = edge->face_right; if ( edge->face_left != floodface ) ids[nids++] = edge->face_left; int result = lwt_be_deleteFacesById(topo, ids, nids); if (result == -1) { _lwt_release_edges(edge, 1); lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } } _lwt_release_edges(edge, 1); return modFace ? floodface : newface.face_id; } LWT_ELEMID lwt_RemEdgeModFace( LWT_TOPOLOGY* topo, LWT_ELEMID edge_id ) { return _lwt_RemEdge( topo, edge_id, 1 ); } LWT_ELEMID lwt_RemEdgeNewFace( LWT_TOPOLOGY* topo, LWT_ELEMID edge_id ) { return _lwt_RemEdge( topo, edge_id, 0 ); } static LWT_ELEMID _lwt_HealEdges( LWT_TOPOLOGY* topo, LWT_ELEMID eid1, LWT_ELEMID eid2, int modEdge ) { LWT_ELEMID ids[2]; LWT_ELEMID commonnode = -1; int caseno = 0; LWT_ISO_EDGE *node_edges; uint64_t num_node_edges; LWT_ISO_EDGE *edges; LWT_ISO_EDGE *e1 = NULL; LWT_ISO_EDGE *e2 = NULL; LWT_ISO_EDGE newedge, updedge, seledge; uint64_t nedges, i; int e1freenode; int e2sign, e2freenode; POINTARRAY *pa; char buf[256]; char *ptr; size_t bufleft = 256; ptr = buf; /* NOT IN THE SPECS: see if the same edge is given twice.. */ if ( eid1 == eid2 ) { lwerror("Cannot heal edge %" LWTFMT_ELEMID " with itself, try with another", eid1); return -1; } ids[0] = eid1; ids[1] = eid2; nedges = 2; edges = lwt_be_getEdgeById(topo, ids, &nedges, LWT_COL_EDGE_ALL); if ((nedges == UINT64_MAX) || (edges == NULL)) { lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } for ( i=0; istart_node == e1->end_node ) { _lwt_release_edges(edges, nedges); lwerror("Edge %" LWTFMT_ELEMID " is closed, cannot heal to edge %" LWTFMT_ELEMID, eid1, eid2); return -1; } if ( e2->start_node == e2->end_node ) { _lwt_release_edges(edges, nedges); lwerror("Edge %" LWTFMT_ELEMID " is closed, cannot heal to edge %" LWTFMT_ELEMID, eid2, eid1); return -1; } /* Find common node */ if ( e1->end_node == e2->start_node ) { commonnode = e1->end_node; caseno = 1; } else if ( e1->end_node == e2->end_node ) { commonnode = e1->end_node; caseno = 2; } /* Check if any other edge is connected to the common node, if found */ if ( commonnode != -1 ) { num_node_edges = 1; node_edges = lwt_be_getEdgeByNode( topo, &commonnode, &num_node_edges, LWT_COL_EDGE_EDGE_ID ); if (num_node_edges == UINT64_MAX) { _lwt_release_edges(edges, nedges); lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } for (i=0; i 0 ) { r = snprintf(ptr, bufleft, "%s%" LWTFMT_ELEMID, ( ptr==buf ? "" : "," ), node_edges[i].edge_id); if ( r >= (int) bufleft ) { bufleft = 0; buf[252] = '.'; buf[253] = '.'; buf[254] = '.'; buf[255] = '\0'; } else { bufleft -= r; ptr += r; } } } lwfree(node_edges); } if ( commonnode == -1 ) { if ( e1->start_node == e2->start_node ) { commonnode = e1->start_node; caseno = 3; } else if ( e1->start_node == e2->end_node ) { commonnode = e1->start_node; caseno = 4; } /* Check if any other edge is connected to the common node, if found */ if ( commonnode != -1 ) { num_node_edges = 1; node_edges = lwt_be_getEdgeByNode( topo, &commonnode, &num_node_edges, LWT_COL_EDGE_EDGE_ID ); if (num_node_edges == UINT64_MAX) { _lwt_release_edges(edges, nedges); lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } for (i=0; i 0 ) { r = snprintf(ptr, bufleft, "%s%" LWTFMT_ELEMID, ( ptr==buf ? "" : "," ), node_edges[i].edge_id); if ( r >= (int) bufleft ) { bufleft = 0; buf[252] = '.'; buf[253] = '.'; buf[254] = '.'; buf[255] = '\0'; } else { bufleft -= r; ptr += r; } } } if ( num_node_edges ) lwfree(node_edges); } } if ( commonnode == -1 ) { _lwt_release_edges(edges, nedges); if ( ptr != buf ) { lwerror("SQL/MM Spatial exception - other edges connected (%s)", buf); } else { lwerror("SQL/MM Spatial exception - non-connected edges"); } return -1; } if ( ! lwt_be_checkTopoGeomRemNode(topo, commonnode, eid1, eid2 ) ) { _lwt_release_edges(edges, nedges); lwerror("%s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } /* Construct the geometry of the new edge */ switch (caseno) { case 1: /* e1.end = e2.start */ pa = ptarray_clone_deep(e1->geom->points); //pa = ptarray_merge(pa, e2->geom->points); ptarray_append_ptarray(pa, e2->geom->points, 0); newedge.start_node = e1->start_node; newedge.end_node = e2->end_node; newedge.next_left = e2->next_left; newedge.next_right = e1->next_right; e1freenode = 1; e2freenode = -1; e2sign = 1; break; case 2: /* e1.end = e2.end */ { POINTARRAY *pa2; pa2 = ptarray_clone_deep(e2->geom->points); ptarray_reverse_in_place(pa2); pa = ptarray_clone_deep(e1->geom->points); //pa = ptarray_merge(e1->geom->points, pa); ptarray_append_ptarray(pa, pa2, 0); ptarray_free(pa2); newedge.start_node = e1->start_node; newedge.end_node = e2->start_node; newedge.next_left = e2->next_right; newedge.next_right = e1->next_right; e1freenode = 1; e2freenode = 1; e2sign = -1; break; } case 3: /* e1.start = e2.start */ pa = ptarray_clone_deep(e2->geom->points); ptarray_reverse_in_place(pa); //pa = ptarray_merge(pa, e1->geom->points); ptarray_append_ptarray(pa, e1->geom->points, 0); newedge.end_node = e1->end_node; newedge.start_node = e2->end_node; newedge.next_left = e1->next_left; newedge.next_right = e2->next_left; e1freenode = -1; e2freenode = -1; e2sign = -1; break; case 4: /* e1.start = e2.end */ pa = ptarray_clone_deep(e2->geom->points); //pa = ptarray_merge(pa, e1->geom->points); ptarray_append_ptarray(pa, e1->geom->points, 0); newedge.end_node = e1->end_node; newedge.start_node = e2->start_node; newedge.next_left = e1->next_left; newedge.next_right = e2->next_right; e1freenode = -1; e2freenode = 1; e2sign = 1; break; default: pa = NULL; e1freenode = 0; e2freenode = 0; e2sign = 0; _lwt_release_edges(edges, nedges); lwerror("Coding error: caseno=%d should never happen", caseno); return -1; break; } newedge.geom = lwline_construct(topo->srid, NULL, pa); if ( modEdge ) { /* Update data of the first edge */ newedge.edge_id = eid1; int result = lwt_be_updateEdgesById(topo, &newedge, 1, LWT_COL_EDGE_NEXT_LEFT | LWT_COL_EDGE_NEXT_RIGHT | LWT_COL_EDGE_START_NODE | LWT_COL_EDGE_END_NODE | LWT_COL_EDGE_GEOM); if (result == -1) { lwline_free(newedge.geom); _lwt_release_edges(edges, nedges); lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } else if (result != 1) { lwline_free(newedge.geom); if ( edges ) _lwt_release_edges(edges, nedges); lwerror("Unexpected error: %d edges updated when expecting 1", i); return -1; } } else { /* Add new edge */ newedge.edge_id = -1; newedge.face_left = e1->face_left; newedge.face_right = e1->face_right; int result = lwt_be_insertEdges(topo, &newedge, 1); if (result == -1) { lwline_free(newedge.geom); _lwt_release_edges(edges, nedges); lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } else if (result == 0) { lwline_free(newedge.geom); _lwt_release_edges(edges, nedges); lwerror("Insertion of split edge failed (no reason)"); return -1; } } lwline_free(newedge.geom); /* -- Update next_left_edge/next_right_edge for -- any edge having them still pointing at the edge being removed -- (eid2 only when modEdge, or both otherwise) -- -- NOTE: -- e#freenode is 1 when edge# end node was the common node -- and -1 otherwise. This gives the sign of possibly found references -- to its "free" (non connected to other edge) endnode. -- e2sign is -1 if edge1 direction is opposite to edge2 direction, -- or 1 otherwise. -- */ /* update edges connected to e2's boundary from their end node */ seledge.next_left = e2freenode * eid2; updedge.next_left = e2freenode * newedge.edge_id * e2sign; int result = lwt_be_updateEdges(topo, &seledge, LWT_COL_EDGE_NEXT_LEFT, &updedge, LWT_COL_EDGE_NEXT_LEFT, NULL, 0); if (result == -1) { _lwt_release_edges(edges, nedges); lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } /* update edges connected to e2's boundary from their start node */ seledge.next_right = e2freenode * eid2; updedge.next_right = e2freenode * newedge.edge_id * e2sign; result = lwt_be_updateEdges(topo, &seledge, LWT_COL_EDGE_NEXT_RIGHT, &updedge, LWT_COL_EDGE_NEXT_RIGHT, NULL, 0); if (result == -1) { _lwt_release_edges(edges, nedges); lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } if ( ! modEdge ) { /* update edges connected to e1's boundary from their end node */ seledge.next_left = e1freenode * eid1; updedge.next_left = e1freenode * newedge.edge_id; result = lwt_be_updateEdges(topo, &seledge, LWT_COL_EDGE_NEXT_LEFT, &updedge, LWT_COL_EDGE_NEXT_LEFT, NULL, 0); if (result == -1) { _lwt_release_edges(edges, nedges); lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } /* update edges connected to e1's boundary from their start node */ seledge.next_right = e1freenode * eid1; updedge.next_right = e1freenode * newedge.edge_id; result = lwt_be_updateEdges(topo, &seledge, LWT_COL_EDGE_NEXT_RIGHT, &updedge, LWT_COL_EDGE_NEXT_RIGHT, NULL, 0); if (result == -1) { _lwt_release_edges(edges, nedges); lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } } /* delete the edges (only second on modEdge or both) */ result = lwt_be_deleteEdges(topo, e2, LWT_COL_EDGE_EDGE_ID); if (result == -1) { _lwt_release_edges(edges, nedges); lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } if ( ! modEdge ) { i = lwt_be_deleteEdges(topo, e1, LWT_COL_EDGE_EDGE_ID); if (result == -1) { _lwt_release_edges(edges, nedges); lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } } _lwt_release_edges(edges, nedges); /* delete the common node */ i = lwt_be_deleteNodesById( topo, &commonnode, 1 ); if (result == -1) { lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } /* -- -- NOT IN THE SPECS: -- Drop composition rows involving second -- edge, as the first edge took its space, -- and all affected TopoGeom have been previously checked -- for being composed by both edges. */ if ( ! lwt_be_updateTopoGeomEdgeHeal(topo, eid1, eid2, newedge.edge_id) ) { lwerror("%s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } return modEdge ? commonnode : newedge.edge_id; } LWT_ELEMID lwt_ModEdgeHeal( LWT_TOPOLOGY* topo, LWT_ELEMID e1, LWT_ELEMID e2 ) { return _lwt_HealEdges( topo, e1, e2, 1 ); } LWT_ELEMID lwt_NewEdgeHeal( LWT_TOPOLOGY* topo, LWT_ELEMID e1, LWT_ELEMID e2 ) { return _lwt_HealEdges( topo, e1, e2, 0 ); } LWT_ELEMID lwt_GetNodeByPoint(LWT_TOPOLOGY *topo, LWPOINT *pt, double tol) { LWT_ISO_NODE *elem; uint64_t num; int flds = LWT_COL_NODE_NODE_ID|LWT_COL_NODE_GEOM; /* geom not needed */ LWT_ELEMID id = 0; POINT2D qp; /* query point */ if ( ! getPoint2d_p(pt->point, 0, &qp) ) { lwerror("Empty query point"); return -1; } elem = lwt_be_getNodeWithinDistance2D(topo, pt, tol, &num, flds, 0); if (num == UINT64_MAX) { lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } else if ( num ) { if ( num > 1 ) { _lwt_release_nodes(elem, num); lwerror("Two or more nodes found"); return -1; } id = elem[0].node_id; _lwt_release_nodes(elem, num); } return id; } LWT_ELEMID lwt_GetEdgeByPoint(LWT_TOPOLOGY *topo, LWPOINT *pt, double tol) { LWT_ISO_EDGE *elem; uint64_t num, i; int flds = LWT_COL_EDGE_EDGE_ID|LWT_COL_EDGE_GEOM; /* GEOM is not needed */ LWT_ELEMID id = 0; LWGEOM *qp = lwpoint_as_lwgeom(pt); /* query point */ if ( lwgeom_is_empty(qp) ) { lwerror("Empty query point"); return -1; } elem = lwt_be_getEdgeWithinDistance2D(topo, pt, tol, &num, flds, 0); if (num == UINT64_MAX) { lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } for (i=0; igeom ) { _lwt_release_edges(elem, num); lwnotice("Corrupted topology: edge %" LWTFMT_ELEMID " has null geometry", e->edge_id); continue; } /* Should we check for intersection not being on an endpoint * as documented ? */ geom = lwline_as_lwgeom(e->geom); dist = lwgeom_mindistance2d_tolerance(geom, qp, tol); if ( dist > tol ) continue; #endif if ( id ) { _lwt_release_edges(elem, num); lwerror("Two or more edges found"); return -1; } else id = e->edge_id; } if ( num ) _lwt_release_edges(elem, num); return id; } LWT_ELEMID lwt_GetFaceByPoint(LWT_TOPOLOGY *topo, const LWPOINT *pt, double tol) { LWT_ELEMID id = 0; LWT_ISO_EDGE *elem; uint64_t num, i; int flds = LWT_COL_EDGE_EDGE_ID | LWT_COL_EDGE_GEOM | LWT_COL_EDGE_FACE_LEFT | LWT_COL_EDGE_FACE_RIGHT; LWGEOM *qp = lwpoint_as_lwgeom(pt); id = lwt_GetFaceContainingPoint(topo, pt); if ( id == -1 ) { lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } if ( id > 0 ) { return id; } if ( tol == 0 ) { return id; } LWDEBUG(1, "No face properly contains query point," " looking for edges"); /* Not in a face, may be in universe or on edge, let's check * for distance */ /* NOTE: we never pass a tolerance of 0 to avoid ever using * ST_Within, which doesn't include endpoints matches */ elem = lwt_be_getEdgeWithinDistance2D(topo, pt, tol?tol:1e-5, &num, flds, 0); if (num == UINT64_MAX) { lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } for (i=0; igeom ) { _lwt_release_edges(elem, num); lwnotice("Corrupted topology: edge %" LWTFMT_ELEMID " has null geometry", e->edge_id); continue; } /* don't consider dangling edges */ if ( e->face_left == e->face_right ) { LWDEBUGF(1, "Edge %" LWTFMT_ELEMID " is dangling, won't consider it", e->edge_id); continue; } geom = lwline_as_lwgeom(e->geom); dist = lwgeom_mindistance2d_tolerance(geom, qp, tol); LWDEBUGF(1, "Distance from edge %" LWTFMT_ELEMID " is %g (tol=%g)", e->edge_id, dist, tol); /* we won't consider edges too far */ if ( dist > tol ) continue; if ( e->face_left == 0 ) { eface = e->face_right; } else if ( e->face_right == 0 ) { eface = e->face_left; } else { _lwt_release_edges(elem, num); lwerror("Two or more faces found"); return -1; } if ( id && id != eface ) { _lwt_release_edges(elem, num); lwerror("Two or more faces found" #if 0 /* debugging */ " (%" LWTFMT_ELEMID " and %" LWTFMT_ELEMID ")", id, eface #endif ); return -1; } else id = eface; } if ( num ) _lwt_release_edges(elem, num); return id; } /* Return the smallest delta that can perturbate * the given value */ static inline double _lwt_minToleranceDouble( double d ) { double ret = 3.6 * pow(10, - ( 15 - log10(d?d:1.0) ) ); return ret; } /* Return the smallest delta that can perturbate * the given point static inline double _lwt_minTolerancePoint2d( const POINT2D* p ) { double max = FP_ABS(p->x); if ( max < FP_ABS(p->y) ) max = FP_ABS(p->y); return _lwt_minToleranceDouble(max); } */ /* Return the smallest delta that can perturbate * the maximum absolute value of a geometry ordinate */ static double _lwt_minTolerance( LWGEOM *g ) { const GBOX* gbox; double max; double ret; gbox = lwgeom_get_bbox(g); if ( ! gbox ) return 0; /* empty */ max = FP_ABS(gbox->xmin); if ( max < FP_ABS(gbox->xmax) ) max = FP_ABS(gbox->xmax); if ( max < FP_ABS(gbox->ymin) ) max = FP_ABS(gbox->ymin); if ( max < FP_ABS(gbox->ymax) ) max = FP_ABS(gbox->ymax); ret = _lwt_minToleranceDouble(max); return ret; } #define _LWT_MINTOLERANCE( topo, geom ) ( \ topo->precision ? topo->precision : _lwt_minTolerance(geom) ) typedef struct scored_pointer_t { void *ptr; double score; } scored_pointer; static int compare_scored_pointer(const void *si1, const void *si2) { double a = ((scored_pointer *)si1)->score; double b = ((scored_pointer *)si2)->score; if ( a < b ) return -1; else if ( a > b ) return 1; else return 0; } /* * @param findFace if non-zero the code will determine which face * contains the given point (unless it is known to be NOT * isolated) * @param moved if not-null will be set to 0 if the point was added * w/out any snapping or 1 otherwise. */ static LWT_ELEMID _lwt_AddPoint(LWT_TOPOLOGY* topo, LWPOINT* point, double tol, int findFace, int *moved) { uint64_t num, i; double mindist = FLT_MAX; LWT_ISO_NODE *nodes, *nodes2; LWT_ISO_EDGE *edges, *edges2; LWGEOM *pt = lwpoint_as_lwgeom(point); int flds; LWT_ELEMID id = 0; scored_pointer *sorted; /* Get tolerance, if 0 was given */ if (!tol) tol = _LWT_MINTOLERANCE(topo, pt); LWDEBUGG(1, pt, "Adding point"); /* -- 1. Check if any existing node is closer than the given precision -- and if so pick the closest TODO: use WithinBox2D */ flds = LWT_COL_NODE_NODE_ID | LWT_COL_NODE_GEOM; nodes = lwt_be_getNodeWithinDistance2D(topo, point, tol, &num, flds, 0); if (num == UINT64_MAX) { lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } if ( num ) { LWDEBUGF(1, "New point is within %.15g units of %d nodes", tol, num); /* Order by distance if there are more than a single return */ if ( num > 1 ) {{ sorted= lwalloc(sizeof(scored_pointer)*num); for (i=0; inode_id, sorted[i].score); } qsort(sorted, num, sizeof(scored_pointer), compare_scored_pointer); nodes2 = lwalloc(sizeof(LWT_ISO_NODE)*num); for (i=0; igeom); double dist = lwgeom_mindistance2d(g, pt); /* TODO: move this check in the previous sort scan ... */ /* must be closer than tolerated, unless distance is zero */ if ( dist && dist >= tol ) continue; if ( ! id || dist < mindist ) { id = n->node_id; mindist = dist; } } if ( id ) { /* found an existing node */ if ( nodes ) _lwt_release_nodes(nodes, num); if ( moved ) *moved = mindist == 0 ? 0 : 1; return id; } } initGEOS(lwnotice, lwgeom_geos_error); /* -- 2. Check if any existing edge falls within tolerance -- and if so split it by a point projected on it TODO: use WithinBox2D */ flds = LWT_COL_EDGE_EDGE_ID|LWT_COL_EDGE_GEOM; edges = lwt_be_getEdgeWithinDistance2D(topo, point, tol, &num, flds, 0); if (num == UINT64_MAX) { lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } if ( num ) { LWDEBUGF(1, "New point is within %.15g units of %d edges", tol, num); /* Order by distance if there are more than a single return */ if ( num > 1 ) {{ int j; sorted = lwalloc(sizeof(scored_pointer)*num); for (i=0; iedge_id, sorted[i].score); } qsort(sorted, num, sizeof(scored_pointer), compare_scored_pointer); edges2 = lwalloc(sizeof(LWT_ISO_EDGE)*num); for (j=0, i=0; igeom); } } num = j; lwfree(sorted); lwfree(edges); edges = edges2; }} for (i=0; igeom); LWGEOM *prj; int contains; LWT_ELEMID edge_id = e->edge_id; LWDEBUGF(1, "Splitting edge %" LWTFMT_ELEMID, edge_id); /* project point to line, split edge by point */ prj = lwgeom_closest_point(g, pt); if ( moved ) *moved = lwgeom_same(prj,pt) ? 0 : 1; if ( lwgeom_has_z(pt) ) {{ /* -- This is a workaround for ClosestPoint lack of Z support: -- http://trac.osgeo.org/postgis/ticket/2033 */ LWGEOM *tmp; double z; POINT4D p4d; LWPOINT *prjpt; /* add Z to "prj" */ tmp = lwgeom_force_3dz(prj, 0); prjpt = lwgeom_as_lwpoint(tmp); getPoint4d_p(point->point, 0, &p4d); z = p4d.z; getPoint4d_p(prjpt->point, 0, &p4d); p4d.z = z; ptarray_set_point4d(prjpt->point, 0, &p4d); lwgeom_free(prj); prj = tmp; }} const POINT2D *pt = getPoint2d_cp(lwgeom_as_lwpoint(prj)->point, 0); contains = ptarray_contains_point_partial(e->geom->points, pt, 0, NULL) == LW_BOUNDARY; if ( ! contains ) {{ double snaptol; LWGEOM *snapedge; LWLINE *snapline; POINT4D p1, p2; LWDEBUGF(1, "Edge %" LWTFMT_ELEMID " does not contain projected point to it", edge_id); /* In order to reduce the robustness issues, we'll pick * an edge that contains the projected point, if possible */ if ( i+1 < num ) { LWDEBUG(1, "But there's another to check"); lwgeom_free(prj); continue; } /* -- The tolerance must be big enough for snapping to happen -- and small enough to snap only to the projected point. -- Unfortunately ST_Distance returns 0 because it also uses -- a projected point internally, so we need another way. */ snaptol = _lwt_minTolerance(prj); snapedge = _lwt_toposnap(g, prj, snaptol); snapline = lwgeom_as_lwline(snapedge); LWDEBUGF(1, "Edge snapped with tolerance %g", snaptol); /* TODO: check if snapping did anything ? */ #if POSTGIS_DEBUG_LEVEL > 0 { size_t sz; char *wkt1 = lwgeom_to_wkt(g, WKT_EXTENDED, 15, &sz); char *wkt2 = lwgeom_to_wkt(snapedge, WKT_EXTENDED, 15, &sz); LWDEBUGF(1, "Edge %s snapped became %s", wkt1, wkt2); lwfree(wkt1); lwfree(wkt2); } #endif /* -- Snapping currently snaps the first point below tolerance -- so may possibly move first point. See ticket #1631 */ getPoint4d_p(e->geom->points, 0, &p1); getPoint4d_p(snapline->points, 0, &p2); LWDEBUGF(1, "Edge first point is %g %g, " "snapline first point is %g %g", p1.x, p1.y, p2.x, p2.y); if ( p1.x != p2.x || p1.y != p2.y ) { LWDEBUG(1, "Snapping moved first point, re-adding it"); if ( LW_SUCCESS != ptarray_insert_point(snapline->points, &p1, 0) ) { lwgeom_free(prj); lwgeom_free(snapedge); _lwt_release_edges(edges, num); lwerror("GEOS exception on Contains: %s", lwgeom_geos_errmsg); return -1; } #if POSTGIS_DEBUG_LEVEL > 0 { size_t sz; char *wkt1 = lwgeom_to_wkt(g, WKT_EXTENDED, 15, &sz); LWDEBUGF(1, "Tweaked snapline became %s", wkt1); lwfree(wkt1); } #endif } #if POSTGIS_DEBUG_LEVEL > 0 else { LWDEBUG(1, "Snapping did not move first point"); } #endif if ( -1 == lwt_ChangeEdgeGeom( topo, edge_id, snapline ) ) { /* TODO: should have invoked lwerror already, leaking memory */ lwgeom_free(prj); lwgeom_free(snapedge); _lwt_release_edges(edges, num); lwerror("lwt_ChangeEdgeGeom failed"); return -1; } lwgeom_free(snapedge); }} #if POSTGIS_DEBUG_LEVEL > 0 else {{ size_t sz; char *wkt1 = lwgeom_to_wkt(g, WKT_EXTENDED, 15, &sz); char *wkt2 = lwgeom_to_wkt(prj, WKT_EXTENDED, 15, &sz); LWDEBUGF(1, "Edge %s contains projected point %s", wkt1, wkt2); lwfree(wkt1); lwfree(wkt2); }} #endif /* TODO: pass 1 as last argument (skipChecks) ? */ id = lwt_ModEdgeSplit( topo, edge_id, lwgeom_as_lwpoint(prj), 0 ); if ( -1 == id ) { /* TODO: should have invoked lwerror already, leaking memory */ lwgeom_free(prj); _lwt_release_edges(edges, num); lwerror("lwt_ModEdgeSplit failed"); return -1; } lwgeom_free(prj); /* * TODO: decimate the two new edges with the given tolerance ? * * the edge identifiers to decimate would be: edge_id and "id" * The problem here is that decimation of existing edges * may introduce intersections or topological inconsistencies, * for example: * * - A node may end up falling on the other side of the edge * - The decimated edge might intersect another existing edge * */ break; /* we only want to snap a single edge */ } _lwt_release_edges(edges, num); } else { /* The point is isolated, add it as such */ /* TODO: pass 1 as last argument (skipChecks) ? */ id = _lwt_AddIsoNode(topo, -1, point, 0, findFace); if ( moved ) *moved = 0; if ( -1 == id ) { /* should have invoked lwerror already, leaking memory */ lwerror("lwt_AddIsoNode failed"); return -1; } } return id; } LWT_ELEMID lwt_AddPoint(LWT_TOPOLOGY* topo, LWPOINT* point, double tol) { return _lwt_AddPoint(topo, point, tol, 1, NULL); } /* Return identifier of an equal edge, 0 if none or -1 on error * (and lwerror gets called on error) * * If an equal edge is found, specify, in "forward" variable whether * the edge is also equal direction-wise * */ static LWT_ELEMID _lwt_GetEqualEdge( LWT_TOPOLOGY *topo, LWLINE *edge, int *forward ) { LWT_ELEMID id; LWT_ISO_EDGE *edges; uint64_t num, i; const GBOX *qbox = lwgeom_get_bbox( lwline_as_lwgeom(edge) ); GEOSGeometry *edgeg; const int flds = LWT_COL_EDGE_EDGE_ID|LWT_COL_EDGE_GEOM; edges = lwt_be_getEdgeWithinBox2D( topo, qbox, &num, flds, 0 ); if (num == UINT64_MAX) { lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } if ( num ) { initGEOS(lwnotice, lwgeom_geos_error); edgeg = LWGEOM2GEOS( lwline_as_lwgeom(edge), 0 ); if ( ! edgeg ) { _lwt_release_edges(edges, num); lwerror("Could not convert edge geometry to GEOS: %s", lwgeom_geos_errmsg); return -1; } for (i=0; igeom); GEOSGeometry *gg; int equals; gg = LWGEOM2GEOS( g, 0 ); if ( ! gg ) { GEOSGeom_destroy(edgeg); _lwt_release_edges(edges, num); lwerror("Could not convert edge geometry to GEOS: %s", lwgeom_geos_errmsg); return -1; } equals = GEOSEquals(gg, edgeg); GEOSGeom_destroy(gg); if ( equals == 2 ) { GEOSGeom_destroy(edgeg); _lwt_release_edges(edges, num); lwerror("GEOSEquals exception: %s", lwgeom_geos_errmsg); return -1; } if ( equals ) { id = e->edge_id; /* Check if direction also matches */ if ( forward ) { /* If input line is closed, we use winding order */ if ( lwline_is_closed(edge) ) { if ( ptarray_isccw(edge->points) == ptarray_isccw(e->geom->points) ) { *forward = 1; } else { *forward = 0; } } else { /* Input line is not closed, checking fist point is enough */ if ( memcmp( getPoint_internal(edge->points, 0), getPoint_internal(e->geom->points, 0), sizeof(POINT2D) ) == 0 ) { *forward = 1; } else { *forward = 0; } } } GEOSGeom_destroy(edgeg); _lwt_release_edges(edges, num); return id; } } GEOSGeom_destroy(edgeg); _lwt_release_edges(edges, num); } return 0; } /* * 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) * * @param forward output parameter, will be populated if * a pre-existing edge was found in the topology, * in which case a value of 1 means the incoming * line will have the same direction of the edge, * and 0 that the incomine line has opposite direction */ static LWT_ELEMID _lwt_AddLineEdge( LWT_TOPOLOGY* topo, LWLINE* edge, double tol, int handleFaceSplit, int *forward ) { LWCOLLECTION *col; LWPOINT *start_point, *end_point; LWGEOM *tmp = 0, *tmp2; LWT_ISO_NODE *node; LWT_ELEMID nid[2]; /* start_node, end_node */ LWT_ELEMID id; /* edge id */ POINT4D p4d; uint64_t nn, i; int moved=0, mm; LWDEBUGG(1, lwline_as_lwgeom(edge), "_lwtAddLineEdge"); LWDEBUGF(1, "_lwtAddLineEdge with tolerance %g", tol); start_point = lwline_get_lwpoint(edge, 0); if ( ! start_point ) { lwnotice("Empty component of noded line"); return 0; /* must be empty */ } nid[0] = _lwt_AddPoint( topo, start_point, _lwt_minTolerance(lwpoint_as_lwgeom(start_point)), handleFaceSplit, &mm ); lwpoint_free(start_point); /* too late if lwt_AddPoint calls lwerror */ if ( nid[0] == -1 ) return -1; /* lwerror should have been called */ moved += mm; end_point = lwline_get_lwpoint(edge, edge->points->npoints-1); if ( ! end_point ) { lwerror("could not get last point of line " "after successfully getting first point !?"); return -1; } nid[1] = _lwt_AddPoint( topo, end_point, _lwt_minTolerance(lwpoint_as_lwgeom(end_point)), handleFaceSplit, &mm ); moved += mm; lwpoint_free(end_point); /* too late if lwt_AddPoint calls lwerror */ if ( nid[1] == -1 ) return -1; /* lwerror should have been called */ /* -- Added endpoints may have drifted due to tolerance, so -- we need to re-snap the edge to the new nodes before adding it */ if ( moved ) { nn = nid[0] == nid[1] ? 1 : 2; node = lwt_be_getNodeById( topo, nid, &nn, LWT_COL_NODE_NODE_ID|LWT_COL_NODE_GEOM ); if (nn == UINT64_MAX) { lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return -1; } start_point = NULL; end_point = NULL; for (i=0; ipoint, 0, &p4d ); lwline_setPoint4d(edge, 0, &p4d); getPoint4d_p( end_point->point, 0, &p4d ); lwline_setPoint4d(edge, edge->points->npoints-1, &p4d); if ( nn ) _lwt_release_nodes(node, nn); /* make valid, after snap (to handle collapses) */ tmp = lwgeom_make_valid(lwline_as_lwgeom(edge)); col = lwgeom_as_lwcollection(tmp); if ( col ) {{ LWCOLLECTION *colex = lwcollection_extract(col, LINETYPE); /* Check if the so-snapped edge collapsed (see #1650) */ if ( colex->ngeoms == 0 ) { lwcollection_free(colex); lwgeom_free(tmp); LWDEBUG(1, "Made-valid snapped edge collapsed"); return 0; } tmp2 = lwgeom_clone_deep(colex->geoms[0]); lwgeom_free(tmp); tmp = tmp2; edge = lwgeom_as_lwline(tmp); lwcollection_free(colex); if ( ! edge ) { /* should never happen */ lwerror("lwcollection_extract(LINETYPE) returned a non-line?"); return -1; } }} else { edge = lwgeom_as_lwline(tmp); if ( ! edge ) { LWDEBUGF(1, "Made-valid snapped edge collapsed to %s", lwtype_name(lwgeom_get_type(tmp))); lwgeom_free(tmp); return 0; } } } /* check if the so-snapped edge _now_ exists */ id = _lwt_GetEqualEdge ( topo, edge, forward ); LWDEBUGF(1, "_lwt_GetEqualEdge returned %" LWTFMT_ELEMID, id); if ( id == -1 ) { if ( tmp ) lwgeom_free(tmp); /* probably too late, due to internal lwerror */ return -1; } if ( id ) { if ( tmp ) lwgeom_free(tmp); /* possibly takes "edge" down with it */ return id; } /* No previously existing edge was found, we'll add one */ /* Remove consecutive vertices below given tolerance * on edge addition */ if ( tol ) {{ tmp2 = lwline_remove_repeated_points(edge, tol); LWDEBUGG(1, tmp2, "Repeated-point removed"); edge = lwgeom_as_lwline(tmp2); if ( tmp ) lwgeom_free(tmp); tmp = tmp2; /* check if the so-decimated edge collapsed to single-point */ if ( nid[0] == nid[1] && edge->points->npoints == 2 ) { lwgeom_free(tmp); LWDEBUG(1, "Repeated-point removed edge collapsed"); return 0; } /* check if the so-decimated edge _now_ exists */ id = _lwt_GetEqualEdge ( topo, edge, forward ); LWDEBUGF(1, "_lwt_GetEqualEdge returned %" LWTFMT_ELEMID, id); if ( id == -1 ) { lwgeom_free(tmp); /* probably too late, due to internal lwerror */ return -1; } if ( id ) { lwgeom_free(tmp); /* takes "edge" down with it */ return id; } }} /* TODO: skip checks ? */ id = _lwt_AddEdge( topo, nid[0], nid[1], edge, 0, handleFaceSplit ? 1 : -1 ); LWDEBUGF(1, "lwt_AddEdgeModFace returned %" LWTFMT_ELEMID, id); if ( id == -1 ) { lwgeom_free(tmp); /* probably too late, due to internal lwerror */ return -1; } lwgeom_free(tmp); /* possibly takes "edge" down with it */ *forward = 1; return id; } /* Simulate split-loop as it was implemented in pl/pgsql version * of TopoGeo_addLinestring */ static LWGEOM * _lwt_split_by_nodes(const LWGEOM *g, const LWGEOM *nodes) { LWCOLLECTION *col = lwgeom_as_lwcollection(nodes); uint32_t i; LWGEOM *bg; bg = lwgeom_clone_deep(g); if ( ! col->ngeoms ) return bg; for (i=0; ingeoms; ++i) { LWGEOM *g2; g2 = lwgeom_split(bg, col->geoms[i]); lwgeom_free(bg); bg = g2; } bg->srid = nodes->srid; return bg; } static LWT_ELEMID* _lwt_AddLine(LWT_TOPOLOGY* topo, LWLINE* line, double tol, int* nedges, int handleFaceSplit) { LWGEOM *geomsbuf[1]; LWGEOM **geoms; uint32_t ngeoms; LWGEOM *noded, *tmp; LWCOLLECTION *col; LWT_ELEMID *ids; LWT_ISO_EDGE *edges; LWT_ISO_NODE *nodes; uint64_t num, numedges = 0, numnodes = 0; uint64_t i; GBOX qbox; int forward; int input_was_closed = 0; POINT4D originalStartPoint; if ( lwline_is_closed(line) ) { input_was_closed = 1; getPoint4d_p( line->points, 0, &originalStartPoint); LWDEBUGF(1, "Input line is closed, original point is %g,%g", originalStartPoint.x, originalStartPoint.y); } *nedges = -1; /* error condition, by default */ /* Get tolerance, if 0 was given */ if ( ! tol ) tol = _LWT_MINTOLERANCE( topo, (LWGEOM*)line ); LWDEBUGF(1, "Working tolerance:%.15g", tol); LWDEBUGF(1, "Input line has srid=%d", line->srid); /* Remove consecutive vertices below given tolerance upfront */ if ( tol ) {{ LWLINE *clean = lwgeom_as_lwline(lwline_remove_repeated_points(line, tol)); tmp = lwline_as_lwgeom(clean); /* NOTE: might collapse to non-simple */ LWDEBUGG(1, tmp, "Repeated-point removed"); }} else tmp=(LWGEOM*)line; /* 1. Self-node */ noded = lwgeom_node((LWGEOM*)tmp); if ( tmp != (LWGEOM*)line ) lwgeom_free(tmp); if ( ! noded ) return NULL; /* should have called lwerror already */ LWDEBUGG(1, noded, "Noded"); qbox = *lwgeom_get_bbox( lwline_as_lwgeom(line) ); LWDEBUGF(1, "Line BOX is %.15g %.15g, %.15g %.15g", qbox.xmin, qbox.ymin, qbox.xmax, qbox.ymax); gbox_expand(&qbox, tol); LWDEBUGF(1, "BOX expanded by %g is %.15g %.15g, %.15g %.15g", tol, qbox.xmin, qbox.ymin, qbox.xmax, qbox.ymax); LWGEOM **nearby = 0; int nearbyindex = 0; int nearbycount = 0; /* 2.0. Find edges falling within tol distance */ edges = lwt_be_getEdgeWithinBox2D( topo, &qbox, &numedges, LWT_COL_EDGE_ALL, 0 ); if (numedges == UINT64_MAX) { lwgeom_free(noded); lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return NULL; } LWDEBUGF(1, "Line has %d points, its bbox intersects %d edges bboxes", line->points->npoints, numedges); if ( numedges ) {{ /* collect those whose distance from us is < tol */ nearbycount += numedges; nearby = lwalloc(numedges * sizeof(LWGEOM *)); for (i=0; igeom); LWDEBUGF(2, "Computing distance from edge %d having %d points", i, e->geom->points->npoints); double dist = lwgeom_mindistance2d(g, noded); /* must be closer than tolerated, unless distance is zero */ if ( dist && dist >= tol ) continue; nearby[nearbyindex++] = g; } LWDEBUGF(1, "Found %d edges closer than tolerance (%g)", nearbyindex, tol); }} int nearbyedgecount = nearbyindex; /* 2.1. Find isolated nodes falling within tol distance * * TODO: add backend-interface support for only getting isolated nodes */ nodes = lwt_be_getNodeWithinBox2D( topo, &qbox, &numnodes, LWT_COL_NODE_ALL, 0 ); if (numnodes == UINT64_MAX) { lwgeom_free(noded); lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return NULL; } LWDEBUGF(1, "Line bbox intersects %d nodes bboxes", numnodes); if ( numnodes ) {{ /* collect those whose distance from us is < tol */ nearbycount = nearbyedgecount + numnodes; nearby = nearby ? lwrealloc(nearby, nearbycount * sizeof(LWGEOM *)) : lwalloc(nearbycount * sizeof(LWGEOM *)) ; int nn = 0; for (i=0; icontaining_face == -1 ) continue; /* skip not-isolated nodes */ LWGEOM *g = lwpoint_as_lwgeom(n->geom); double dist = lwgeom_mindistance2d(g, noded); /* must be closer than tolerated, unless distance is zero */ if ( dist && dist >= tol ) { LWDEBUGF(1, "Node %d is %g units away, we tolerate only %g", n->node_id, dist, tol); continue; } nearby[nearbyindex++] = g; nn = nn + 1; } LWDEBUGF(1, "Found %d isolated nodes closer than tolerance (%g)", nn, tol); }} int nearbynodecount = nearbyindex - nearbyedgecount; nearbycount = nearbyindex; LWDEBUGF(1, "Number of nearby elements is %d", nearbycount); /* 2.2. Snap to nearby elements */ if ( nearbycount ) {{ LWCOLLECTION *col; LWGEOM *elems; col = lwcollection_construct(COLLECTIONTYPE, topo->srid, NULL, nearbycount, nearby); elems = lwcollection_as_lwgeom(col); LWDEBUGG(1, elems, "Collected nearby elements"); tmp = _lwt_toposnap(noded, elems, tol); lwgeom_free(noded); noded = tmp; LWDEBUGG(1, noded, "Elements-snapped"); if ( input_was_closed ) {{ /* Recompute start point in case it moved */ LWLINE *scrolled = lwgeom_as_lwline(noded); if (scrolled) { getPoint4d_p( scrolled->points, 0, &originalStartPoint); LWDEBUGF(1, "Closed input line start point after snap %g,%g", originalStartPoint.x, originalStartPoint.y); } }} /* will not release the geoms array */ lwcollection_release(col); /* -- re-node to account for ST_Snap introduced self-intersections -- See http://trac.osgeo.org/postgis/ticket/1714 -- TODO: consider running UnaryUnion once after all noding */ tmp = lwgeom_unaryunion(noded); lwgeom_free(noded); noded = tmp; LWDEBUGG(1, noded, "Unary-unioned"); }} /* 2.3. Node with nearby edges */ if ( nearbyedgecount ) {{ LWCOLLECTION *col; LWGEOM *iedges; /* just an alias for col */ LWGEOM *diff, *xset; LWDEBUGF(1, "Line intersects %d edges", nearbyedgecount); col = lwcollection_construct(COLLECTIONTYPE, topo->srid, NULL, nearbyedgecount, nearby); iedges = lwcollection_as_lwgeom(col); LWDEBUGG(1, iedges, "Collected edges"); LWDEBUGF(1, "Diffing noded, with srid=%d " "and interesecting edges, with srid=%d", noded->srid, iedges->srid); diff = lwgeom_difference(noded, iedges); LWDEBUGG(1, diff, "Differenced"); LWDEBUGF(1, "Intersecting noded, with srid=%d " "and interesecting edges, with srid=%d", noded->srid, iedges->srid); xset = lwgeom_intersection(noded, iedges); LWDEBUGG(1, xset, "Intersected"); lwgeom_free(noded); /* We linemerge here because INTERSECTION, as of GEOS 3.8, * will result in shared segments being output as multiple * lines rather than a single line. Example: INTERSECTION( 'LINESTRING(0 0, 5 0, 8 0, 10 0,12 0)', 'LINESTRING(5 0, 8 0, 10 0)' ) == MULTILINESTRING((5 0,8 0),(8 0,10 0)) * We will re-split in a subsequent step, by splitting * the final line with pre-existing nodes */ LWDEBUG(1, "Linemerging intersection"); tmp = lwgeom_linemerge(xset); LWDEBUGG(1, tmp, "Linemerged"); lwgeom_free(xset); xset = tmp; /* * Here we union the (linemerged) intersection with * the difference (new lines) */ LWDEBUG(1, "Unioning difference and (linemerged) intersection"); noded = lwgeom_union(diff, xset); LWDEBUGG(1, noded, "Diff-Xset Unioned"); lwgeom_free(xset); lwgeom_free(diff); /* will not release the geoms array */ lwcollection_release(col); if ( input_was_closed ) {{ LWLINE *scrolled = lwgeom_as_lwline(noded); if (scrolled) { if ( lwline_is_closed(scrolled) ) { ptarray_scroll_in_place(scrolled->points, &originalStartPoint); } else { LWDEBUGG(1, lwline_as_lwgeom(scrolled), "Linemerged intersected input is not closed anymore"); } } else { LWDEBUGG(1, xset, "Linemerged intersected input is not a line anymore"); } }} }} /* 2.4. Split by pre-existing nodes * * Pre-existing nodes are isolated nodes AND endpoints * of intersecting edges */ if ( nearbyedgecount ) { nearbycount += nearbyedgecount * 2; /* make space for endpoints */ nearby = lwrealloc(nearby, nearbycount * sizeof(LWGEOM *)); for (int i=0; ipoints->npoints-1); /* TODO: only add if within distance to noded AND if not duplicated */ nearby[nearbyindex++] = lwpoint_as_lwgeom(startNode); nearbynodecount++; nearby[nearbyindex++] = lwpoint_as_lwgeom(endNode); nearbynodecount++; } } if ( nearbynodecount ) { col = lwcollection_construct(MULTIPOINTTYPE, topo->srid, NULL, nearbynodecount, nearby + nearbyedgecount); LWGEOM *inodes = lwcollection_as_lwgeom(col); /* TODO: use lwgeom_split of lwgeom_union ... */ tmp = _lwt_split_by_nodes(noded, inodes); lwgeom_free(noded); noded = tmp; LWDEBUGG(1, noded, "Node-split"); /* will not release the geoms array */ lwcollection_release(col); } LWDEBUG(1, "Freeing up nearby elements"); /* TODO: free up endpoints of nearbyedges */ if ( nearby ) lwfree(nearby); if ( nodes ) _lwt_release_nodes(nodes, numnodes); if ( edges ) _lwt_release_edges(edges, numedges); LWDEBUGG(1, noded, "Finally-noded"); /* 3. For each (now-noded) segment, insert an edge */ col = lwgeom_as_lwcollection(noded); if ( col ) { LWDEBUG(1, "Noded line was a collection"); geoms = col->geoms; ngeoms = col->ngeoms; } else { LWDEBUG(1, "Noded line was a single geom"); geomsbuf[0] = noded; geoms = geomsbuf; ngeoms = 1; } LWDEBUGF(1, "Line was split into %d edges", ngeoms); /* TODO: refactor to first add all nodes (re-snapping edges if * needed) and then check all edges for existing already * ( so to save a DB scan for each edge to be added ) */ ids = lwalloc(sizeof(LWT_ELEMID)*ngeoms); num = 0; for ( i=0; isrid = noded->srid; #if POSTGIS_DEBUG_LEVEL > 0 { size_t sz; char *wkt1 = lwgeom_to_wkt(g, WKT_EXTENDED, 15, &sz); LWDEBUGF(1, "Component %d of split line is: %s", i, wkt1); lwfree(wkt1); } #endif id = _lwt_AddLineEdge( topo, lwgeom_as_lwline(g), tol, handleFaceSplit, &forward ); LWDEBUGF(1, "_lwt_AddLineEdge returned %" LWTFMT_ELEMID, id); if ( id < 0 ) { lwgeom_free(noded); lwfree(ids); return NULL; } if ( ! id ) { LWDEBUGF(1, "Component %d of split line collapsed", i); continue; } LWDEBUGF(1, "Component %d of split line is %s edge %" LWTFMT_ELEMID, i, forward ? "forward" : "backward", id); ids[num++] = forward ? id : -id; /* TODO: skip duplicates */ } LWDEBUGG(1, noded, "Noded before free"); lwgeom_free(noded); /* TODO: XXX remove duplicated ids if not done before */ *nedges = num; 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) { uint32_t i; *nfaces = -1; /* error condition, by default */ int num; LWT_ISO_FACE *faces; uint64_t nfacesinbox; uint64_t j; LWT_ELEMID *ids = NULL; GBOX qbox; const GEOSPreparedGeometry *ppoly; GEOSGeometry *polyg; /* Get tolerance, if 0 was given */ if ( ! tol ) tol = _LWT_MINTOLERANCE( topo, (LWGEOM*)poly ); LWDEBUGF(1, "Working tolerance:%.15g", tol); /* Add each ring as an edge */ for ( i=0; inrings; ++i ) { LWLINE *line; POINTARRAY *pa; LWT_ELEMID *eids; int nedges; pa = ptarray_clone(poly->rings[i]); line = lwline_construct(topo->srid, NULL, pa); eids = lwt_AddLine( topo, line, tol, &nedges ); if ( nedges < 0 ) { /* probably too late as lwt_AddLine invoked lwerror */ lwline_free(line); lwerror("Error adding ring %d of polygon", i); return NULL; } lwline_free(line); lwfree(eids); } /* -- Find faces covered by input polygon -- NOTE: potential snapping changed polygon edges */ qbox = *lwgeom_get_bbox( lwpoly_as_lwgeom(poly) ); gbox_expand(&qbox, tol); faces = lwt_be_getFaceWithinBox2D( topo, &qbox, &nfacesinbox, LWT_COL_FACE_ALL, 0 ); if (nfacesinbox == UINT64_MAX) { lwfree(ids); lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); return NULL; } num = 0; if ( nfacesinbox ) { polyg = LWGEOM2GEOS(lwpoly_as_lwgeom(poly), 0); if ( ! polyg ) { _lwt_release_faces(faces, nfacesinbox); lwerror("Could not convert poly geometry to GEOS: %s", lwgeom_geos_errmsg); return NULL; } ppoly = GEOSPrepare(polyg); ids = lwalloc(sizeof(LWT_ELEMID)*nfacesinbox); for ( j=0; jface_id ); if ( ! fg ) { j = f->face_id; /* so we can destroy faces */ GEOSPreparedGeom_destroy(ppoly); GEOSGeom_destroy(polyg); lwfree(ids); _lwt_release_faces(faces, nfacesinbox); lwerror("Could not get geometry of face %" LWTFMT_ELEMID, j); return NULL; } /* check if a point on this face's surface is covered by our polygon */ fgg = LWGEOM2GEOS(fg, 0); lwgeom_free(fg); if ( ! fgg ) { GEOSPreparedGeom_destroy(ppoly); GEOSGeom_destroy(polyg); _lwt_release_faces(faces, nfacesinbox); lwerror("Could not convert edge geometry to GEOS: %s", lwgeom_geos_errmsg); return NULL; } sp = GEOSPointOnSurface(fgg); GEOSGeom_destroy(fgg); if ( ! sp ) { GEOSPreparedGeom_destroy(ppoly); GEOSGeom_destroy(polyg); _lwt_release_faces(faces, nfacesinbox); lwerror("Could not find point on face surface: %s", lwgeom_geos_errmsg); return NULL; } covers = GEOSPreparedCovers( ppoly, sp ); GEOSGeom_destroy(sp); if (covers == 2) { GEOSPreparedGeom_destroy(ppoly); GEOSGeom_destroy(polyg); _lwt_release_faces(faces, nfacesinbox); lwerror("PreparedCovers error: %s", lwgeom_geos_errmsg); return NULL; } if ( ! covers ) { continue; /* we're not composed by this face */ } /* TODO: avoid duplicates ? */ ids[num++] = f->face_id; } GEOSPreparedGeom_destroy(ppoly); GEOSGeom_destroy(polyg); _lwt_release_faces(faces, nfacesinbox); } /* possibly 0 if non face's surface point was found * to be covered by input polygon */ *nfaces = num; 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 >= (int) 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 == NULL || 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(__attribute__((__unused__)) 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; uint64_t nelems = 1; edge = lwt_be_getEdgeWithinBox2D( topo, NULL, &nelems, fields, 0); *numedges = nelems; if (nelems == UINT64_MAX) { 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(__attribute__((__unused__)) 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; /* 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)", edge->edge_id * side, 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)", edge->edge_id * side, 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; uint64_t 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 == UINT64_MAX) { 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; ibe_iface)); /* cberror(topo->be_data, "Error from cb_getClosestEdge"); */ return -1; } if (numedges == 0) { /* If there are no edges the point is in the universal face */ return 0; } LWDEBUGGF(2, lwline_as_lwgeom(closestEdge->geom), "Closest edge %" LWTFMT_ELEMID, closestEdge->edge_id); /* Find closest segment of edge to the point */ queryPoint = getPoint2d_cp(pt->point, 0); closestSegmentIndex = ptarray_closest_segment_2d(closestEdge->geom->points, queryPoint, &dist); LWDEBUGF(1, "Closest segment on edge %" LWTFMT_ELEMID " is %d (dist %g)", closestEdge->edge_id, closestSegmentIndex, dist); closestSegmentP0 = getPoint2d_cp(closestEdge->geom->points, closestSegmentIndex); closestSegmentP1 = getPoint2d_cp(closestEdge->geom->points, closestSegmentIndex + 1); LWDEBUGF(1, "Closest segment on edge %" LWTFMT_ELEMID " is LINESTRING(%g %g, %g %g)", closestEdge->edge_id, closestSegmentP0->x, closestSegmentP0->y, closestSegmentP1->x, closestSegmentP1->y ); /* * We use comp.graphics.algorithms Frequently Asked Questions method * * (1) AC dot AB * r = ---------- * ||AB||^2 * r has the following meaning: * r=0 P = A * r=1 P = B * r<0 P is on the backward extension of AB * r>1 P is on the forward extension of AB * 0x-A->x) * (B->x-A->x) + (p->y-A->y) * (B->y-A->y) )/( (B->x-A->x)*(B->x-A->x) +(B->y-A->y)*(B->y-A->y) ); if ( r <= 0 ) { closestPointOnEdge = A; closestPointVertex = closestSegmentIndex; if ( closestSegmentIndex == 0 ) { closestNode = closestEdge->start_node; } } else if (r >= 1 ) { closestPointOnEdge = B; closestPointVertex = closestSegmentIndex + 1; if ( closestSegmentIndex + 2 == closestEdge->geom->points->npoints ) { closestNode = closestEdge->end_node; } } else { closestPointVertex = closestEdge->geom->points->npoints; } if ( closestNode != 0 ) { LWDEBUGF(1, "Closest point is node %d", closestNode); if ( dist == 0 ) { LWDEBUGF(1, "Query point is node %d", closestNode); /* Query point is the node * * If all edges incident to the node are * dangling, we can return their common * side face, otherwise the point will be * on multiple face boundaries */ if ( closestEdge->face_left != closestEdge->face_right ) { _lwt_release_edges(closestEdge, 1); lwerror("Two or more faces found"); return -1; } containingFace = closestEdge->face_left; /* Check other incident edges */ numedges = 1; edges = lwt_be_getEdgeByNode( topo, &closestNode, &numedges, LWT_COL_EDGE_FACE_LEFT|LWT_COL_EDGE_FACE_RIGHT ); if (numedges == UINT64_MAX) { lwerror("Backend error from getEdgeByNode: %s", lwt_be_lastErrorMessage(topo->be_iface)); /* cberror(topo->be_data, "Error from cb_getClosestEdge"); */ _lwt_release_edges(closestEdge, 1); return -1; } for (i=0; iedge_id); _lwt_release_edges(edges, numedges); _lwt_release_edges(closestEdge, 1); return -1; } LWDEBUGF(1, "lwt_be_getEdgeByNode returned %d edges", numedges); _lwt_release_edges(edges, numedges); _lwt_release_edges(closestEdge, 1); return containingFace; } /* Closest point is a node, but query point is NOT on the node */ /* let's do azimuth computation */ edgeend ee; if ( ! azimuth_pt_pt(closestPointOnEdge, queryPoint, &ee.myaz) ) { lwerror("error computing azimuth of query point [%.15g %.15g,%.15g %.15g]", closestPointOnEdge->x, closestPointOnEdge->y, queryPoint->x, queryPoint->y); _lwt_release_edges(closestEdge, 1); return -1; } LWDEBUGF(1, "Query point azimuth is %g", ee.myaz); int found = _lwt_FindAdjacentEdges( topo, closestNode, &ee, NULL, -1 ); if ( ! found ) { lwerror("Unexpected backend return: _lwt_FindAdjacentEdges(%d) found no edges when we previously found edge %d ending on that node", closestNode, closestEdge->edge_id); _lwt_release_edges(closestEdge, 1); return -1; } _lwt_release_edges(closestEdge, 1); return ee.cwFace; } LWDEBUG(1, "Closest point is NOT a node"); /* If this edge has the same face on the left and right sides * we found the face containing our query point */ if ( closestEdge->face_left == closestEdge->face_right ) { containingFace = closestEdge->face_left; _lwt_release_edges(closestEdge, 1); return containingFace; } if ( dist == 0 ) { /* We checked the dangling case above */ _lwt_release_edges(closestEdge, 1); lwerror("Two or more faces found"); return -1; } /* Find on which side of the segment the query point lays */ if ( closestPointVertex != closestEdge->geom->points->npoints ) { /* Closest point is a vertex of the closest segment */ LWDEBUGF(1, "Closest point is vertex %d of closest segment", closestPointVertex); /* * We need to check if rotating clockwise the line * from previous vertex to closest vertex clockwise * around the closest vertex encounters * the line to query point first (which means it's on the left * of the closest edge) or the line to next vertex first (which * means the query point is on the right) */ uint32_t prevVertexIndex = closestPointVertex > 0 ? closestPointVertex - 1u : closestEdge->geom->points->npoints - 2u; /* last vertex would be == first vertex, this being a closed edge */ const POINT2D *prevVertex = getPoint2d_cp( closestEdge->geom->points, prevVertexIndex ); LWDEBUGF(1, "Previous vertex %u is POINT(%.15g %.15g)", prevVertexIndex, prevVertex->x, prevVertex->y ); uint32_t nextVertexIndex = closestPointVertex == closestEdge->geom->points->npoints - 1u ? 1u : /* first point would be == last point, this being a closed edge */ closestPointVertex + 1u; const POINT2D *nextVertex = getPoint2d_cp( closestEdge->geom->points, nextVertexIndex ); LWDEBUGF(1, "Next vertex %u is POINT(%.15g %.15g)", nextVertexIndex, nextVertex->x, nextVertex->y ); double azS0; /* azimuth from previous vertex to closestPointVertex */ double azS1; /* azimuth from closestPointVertex to next vertex */ double azSL; /* azimuth from closestPointVertex to query point */ if ( ! azimuth_pt_pt(closestPointOnEdge, prevVertex, &azS0)) { lwerror("error computing azimuth of segment to closest point [%.15g %.15g,%.15g %.15g]", closestPointOnEdge->x, closestPointOnEdge->y, prevVertex->x, prevVertex->y); _lwt_release_edges(closestEdge, 1); return -1; } if ( ! azimuth_pt_pt(closestPointOnEdge, nextVertex, &azS1)) { lwerror("error computing azimuth of segment from closest point [%.15g %.15g,%.15g %.15g]", closestPointOnEdge->x, closestPointOnEdge->y, nextVertex->x, nextVertex->y); _lwt_release_edges(closestEdge, 1); return -1; } if ( ! azimuth_pt_pt(closestPointOnEdge, queryPoint, &azSL) ) { lwerror("error computing azimuth of queryPoint [%.15g %.15g,%.15g %.15g]", closestPointOnEdge->x, closestPointOnEdge->y, queryPoint->x, queryPoint->y); _lwt_release_edges(closestEdge, 1); return -1; } double angle_S0_S1 = azS1 - azS0; if ( angle_S0_S1 < 0 ) angle_S0_S1 += 2 * M_PI; double angle_S0_SL = azSL - azS0; if ( angle_S0_SL < 0 ) angle_S0_SL += 2 * M_PI; LWDEBUGF(1, "Azimuths from closest (vertex) point: P:%g, N:%g (+%g), Q:%g (+%g)", azS0, azS1, angle_S0_S1, azSL, angle_S0_SL ); if ( angle_S0_SL < angle_S0_S1 ) { /* line to query point was encountered first, is on the left */ containingFace = closestEdge->face_left; } else { /* line to query point was NOT encountered first, is on the right */ containingFace = closestEdge->face_right; } } else { /* Closest point is internal to closest segment, we can use * lw_segment_side */ LWDEBUGF(1, "Closest point is internal to closest segment, calling lw_segment_side((%g,%g),(%g,%g),(%g,%g)", closestSegmentP0->x, closestSegmentP0->y, closestSegmentP1->x, closestSegmentP1->y, queryPoint->x, queryPoint->y ); closestSegmentSide = lw_segment_side(closestSegmentP0, closestSegmentP1, queryPoint); LWDEBUGF(1, "Side of closest segment query point falls on: %d", closestSegmentSide); if ( closestSegmentSide == -1 ) /* left */ { containingFace = closestEdge->face_left; } else if ( closestSegmentSide == 1 ) /* right */ { containingFace = closestEdge->face_right; } else { lwerror("Unexpected collinearity reported from lw_segment_side"); _lwt_release_edges(closestEdge, 1); return -1; } } _lwt_release_edges(closestEdge, 1); return containingFace; }