From: Sandro Santilli Date: Thu, 20 Aug 2015 10:30:13 +0000 (+0000) Subject: Implement TopoGeo_AddLinestring to C X-Git-Tag: 2.2.0rc1~140 X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=bc6909f2d23c6ac742b1e6b91e2504e45fb3a33a;p=postgis Implement TopoGeo_AddLinestring to C Also: - Convert srid=-1 in topology to officially unknown srid at load time (See #3192). - Use hexwkb for box-based callback queries to avoid drifts. - Fix minTolerance computation. Funded by Tuscany Region (Italy) - SITA (CIG: 60351023B8) git-svn-id: http://svn.osgeo.org/postgis/trunk@13948 b70326c6-7e19-0410-871a-916f4a2858ee --- diff --git a/liblwgeom/TODO b/liblwgeom/TODO index ceca0b24c..2c6f4ce31 100644 --- a/liblwgeom/TODO +++ b/liblwgeom/TODO @@ -26,5 +26,5 @@ lwt_GetFaceByPoint X lwt_GetNodeByPoint X lwt_GetEdgeByPoint X lwt_TopoGeo_AddPoint X -lwt_TopoGeo_AddLineString +lwt_TopoGeo_AddLineString X lwt_TopoGeo_AddPolygon diff --git a/liblwgeom/lwgeom_topo.c b/liblwgeom/lwgeom_topo.c index 673c6ea19..6636a8f56 100644 --- a/liblwgeom/lwgeom_topo.c +++ b/liblwgeom/lwgeom_topo.c @@ -18,7 +18,7 @@ #include "../postgis_config.h" -/*#define POSTGIS_DEBUG_LEVEL 1*/ +#define POSTGIS_DEBUG_LEVEL 1 #include "lwgeom_log.h" #include "liblwgeom_internal.h" @@ -36,6 +36,17 @@ # define LWTFMT_ELEMID PRId64 #endif +/* TODO: move this to lwgeom_log.h */ +#define LWDEBUGG(level, geom, msg, ...) \ + if (POSTGIS_DEBUG_LEVEL >= level) \ + do { \ + size_t sz; \ + char *wkt1 = lwgeom_to_wkt(geom, WKT_ISO, 15, &sz); \ + LWDEBUGF(level, msg ": %s", __VA_ARGS__ wkt1); \ + lwfree(wkt1); \ + } while (0); + + /********************************************************************* * * Backend iface @@ -4593,6 +4604,7 @@ _lwt_minTolerance( LWGEOM *g ) { const GBOX* gbox; double max; + double ret; gbox = lwgeom_get_bbox(g); if ( ! gbox ) return 0; /* empty */ @@ -4601,7 +4613,9 @@ _lwt_minTolerance( LWGEOM *g ) if ( max < FP_ABS(gbox->ymin) ) max = FP_ABS(gbox->ymin); if ( max < FP_ABS(gbox->ymax) ) max = FP_ABS(gbox->ymax); - return 3.6 * pow(10, - ( 15 - log(max?max:1.0) ) ); + ret = 3.6 * pow(10, - ( 15 - log10(max?max:1.0) ) ); + + return ret; } #define _LWT_MINTOLERANCE( topo, geom ) ( \ @@ -4618,6 +4632,7 @@ lwt_AddPoint(LWT_TOPOLOGY* topo, LWPOINT* point, double tol) int flds; LWT_ELEMID id = 0; + /* Get tolerance, if 0 was given */ if ( ! tol ) tol = _LWT_MINTOLERANCE( topo, pt ); /* @@ -4849,12 +4864,488 @@ lwt_AddPoint(LWT_TOPOLOGY* topo, LWPOINT* point, double tol) return id; } +/* Return identifier of an equal edge, 0 if none or -1 on error + * (and lwerror gets called on error) + */ +static LWT_ELEMID +_lwt_GetEqualEdge( LWT_TOPOLOGY *topo, LWLINE *edge ) +{ + LWT_ELEMID id; + LWT_ISO_EDGE *edges; + int 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 == -1 ) + { + 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; + 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-splitted line edge. Used by lwt_AddLine + * Return edge id, 0 if none added (empty edge), -1 on error + */ +static LWT_ELEMID +_lwt_AddLineEdge( LWT_TOPOLOGY* topo, LWLINE* edge, double tol ) +{ + LWCOLLECTION *col; + LWPOINT *start_point, *end_point; + LWGEOM *tmp; + LWT_ISO_NODE *node; + LWT_ELEMID nid[2]; /* start_node, end_node */ + LWT_ELEMID id; /* edge id */ + POINT4D p4d; + int nn, i; + + 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, tol ); + lwpoint_free(start_point); /* too late if lwt_AddPoint calls lwerror */ + if ( nid[0] == -1 ) return -1; /* lwerror should have been called */ + + 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, tol ); + 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 + */ + + 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 == -1 ) + { + 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 ) + {{ + LWGEOM *tmp2; + + col = lwcollection_extract(col, LINETYPE); + + /* Check if the so-snapped edge collapsed (see #1650) */ + if ( col->ngeoms == 0 ) + { + lwcollection_free(col); + lwgeom_free(tmp); + LWDEBUG(1, "Made-valid snapped edge collapsed"); + return 0; + } + + tmp2 = lwgeom_clone_deep( col->geoms[0] ); + lwgeom_free(tmp); + tmp = tmp2; + edge = lwgeom_as_lwline(tmp); + lwcollection_free(col); + 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 ); + 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); /* possibly takes "edge" down with it */ + return id; + } + + /* No previously existing edge was found, we'll add one */ + /* TODO: skip checks, I guess ? */ + id = lwt_AddEdgeModFace( topo, nid[0], nid[1], edge, 0 ); + 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 */ + + 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); + int 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; + } + + return bg; +} + LWT_ELEMID* lwt_AddLine(LWT_TOPOLOGY* topo, LWLINE* line, double tol, int* nedges) { - *nedges = -1; - lwerror("Not implemented yet"); - return NULL; + LWGEOM *geomsbuf[1]; + LWGEOM **geoms; + int ngeoms; + LWGEOM *noded; + LWCOLLECTION *col; + LWT_ELEMID *ids; + LWT_ISO_EDGE *edges; + LWT_ISO_NODE *nodes; + int num; + int i; + GBOX qbox; + + *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); + + /* 1. Self-node */ + noded = lwgeom_node((LWGEOM*)line); + 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); + + /* 2. Node to edges falling within tol distance */ + edges = lwt_be_getEdgeWithinBox2D( topo, &qbox, &num, LWT_COL_EDGE_ALL, 0 ); + if ( num == -1 ) + { + lwgeom_free(noded); + lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); + return NULL; + } + LWDEBUGF(1, "Line bbox intersects %d edges bboxes", num); + if ( num ) + {{ + /* collect those whose distance from us is < tol */ + LWGEOM **nearby = lwalloc(sizeof(LWGEOM *)*num); + int nn=0; + for (i=0; igeom); + double dist = lwgeom_mindistance2d(g, noded); + if ( dist >= tol ) continue; /* must be closer than tolerated */ + nearby[nn++] = g; + } + if ( nn ) + {{ + LWCOLLECTION *col; + LWGEOM *iedges; /* just an alias for col */ + LWGEOM *snapped; + LWGEOM *set1, *set2; + + LWDEBUGF(1, "Line intersects %d edges", nn); + + col = lwcollection_construct(COLLECTIONTYPE, topo->srid, + NULL, nn, nearby); + iedges = lwcollection_as_lwgeom(col); + LWDEBUGG(1, iedges, "Collected edges"); + LWDEBUGF(1, "Snapping noded, with srid=%d " + "to interesecting edges, with srid=%d", + noded->srid, iedges->srid); + snapped = lwgeom_snap(noded, iedges, tol); + lwgeom_free(noded); + LWDEBUGG(1, snapped, "Snapped"); + LWDEBUGF(1, "Diffing snapped, with srid=%d " + "and interesecting edges, with srid=%d", + snapped->srid, iedges->srid); + noded = lwgeom_difference(snapped, iedges); + LWDEBUGG(1, noded, "Differenced"); + LWDEBUGF(1, "Intersecting snapped, with srid=%d " + "and interesecting edges, with srid=%d", + snapped->srid, iedges->srid); + set1 = lwgeom_intersection(snapped, iedges); + LWDEBUGG(1, set1, "Intersected"); + lwgeom_free(snapped); + LWDEBUGF(1, "Linemerging set1, with srid=%d", set1->srid); + set2 = lwgeom_linemerge(set1); + LWDEBUGG(1, set2, "Linemerged"); + LWDEBUGG(1, noded, "Noded"); + lwgeom_free(set1); + LWDEBUGF(1, "Unioning noded, with srid=%d " + "and set2, with srid=%d", noded->srid, set2->srid); + set1 = lwgeom_union(noded, set2); + lwgeom_free(set2); + lwgeom_free(noded); + noded = set1; + LWDEBUGG(1, set1, "Unioned"); + + /* will not release the geoms array */ + lwcollection_release(col); + }} + lwfree(nearby); + _lwt_release_edges(edges, num); + }} + + /* 2.1. Node with existing nodes within tol + * TODO: check if we should be only considering _isolated_ nodes! */ + nodes = lwt_be_getNodeWithinBox2D( topo, &qbox, &num, LWT_COL_EDGE_ALL, 0 ); + if ( num == -1 ) + { + lwgeom_free(noded); + lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface)); + return NULL; + } + LWDEBUGF(1, "Line bbox intersects %d nodes bboxes", num); + if ( num ) + {{ + /* collect those whose distance from us is < tol */ + LWGEOM **nearby = lwalloc(sizeof(LWGEOM *)*num); + int nn=0; + for (i=0; igeom); + double dist = lwgeom_mindistance2d(g, noded); + if ( dist >= tol ) continue; /* must be closer than tolerated */ + nearby[nn++] = g; + } + if ( nn ) + {{ + LWCOLLECTION *col; + LWGEOM *inodes; /* just an alias for col */ + LWGEOM *tmp; + + LWDEBUGF(1, "Line intersects %d nodes", nn); + + col = lwcollection_construct(MULTIPOINTTYPE, topo->srid, + NULL, nn, nearby); + inodes = lwcollection_as_lwgeom(col); + + LWDEBUGG(1, inodes, "Collected nodes"); + + /* TODO: consider snapping once against all elements + * (rather than once with edges and once with nodes) */ + tmp = lwgeom_snap(noded, inodes, tol); + lwgeom_free(noded); + noded = tmp; + LWDEBUGG(1, noded, "Node-snapped"); + + tmp = _lwt_split_by_nodes(noded, inodes); + /* lwgeom_split(noded, inodes); */ + lwgeom_free(noded); + noded = tmp; + LWDEBUGG(1, noded, "Node-splitted"); + + /* 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"); + + }} + lwfree(nearby); + _lwt_release_nodes(nodes, num); + }} + + LWDEBUG(1, "XXX"); + { size_t sz; + lwnotice("%s", lwgeom_to_wkt(noded, WKT_ISO, 15, &sz)); } + 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 splitted 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; i 0 + { + size_t sz; + char *wkt1 = lwgeom_to_wkt(g, WKT_ISO, 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 ); + 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 edge %" LWTFMT_ELEMID, + i, id); + ids[num++] = id; /* TODO: skip duplicates */ + } + + LWDEBUG(1, "YYY"); + { size_t sz; + lwnotice("%s", lwgeom_to_wkt(noded, WKT_ISO, 15, &sz)); } + LWDEBUGG(1, noded, "Noded before free"); + lwgeom_release(noded); /* for some reason freeing this here errors out */ + + /* TODO: XXX remove duplicated ids if not done before */ + + *nedges = num; + return ids; } LWT_ELEMID* diff --git a/topology/postgis_topology.c b/topology/postgis_topology.c index ae55fe01a..b71525154 100644 --- a/topology/postgis_topology.c +++ b/topology/postgis_topology.c @@ -24,7 +24,7 @@ #include "liblwgeom_internal.h" /* for gbox_clone */ #include "liblwgeom_topo.h" -/*#define POSTGIS_DEBUG_LEVEL 1*/ +#define POSTGIS_DEBUG_LEVEL 1 #include "lwgeom_log.h" #include "lwgeom_pg.h" @@ -94,6 +94,41 @@ cberror(const LWT_BE_DATA* be_in, const char *fmt, ...) va_end(ap); } +static void +_lwtype_upper_name(int type, char *buf, size_t buflen) +{ + char *ptr; + snprintf(buf, buflen, "%s", lwtype_name(type)); + buf[buflen-1] = '\0'; + ptr = buf; + while (*ptr) { + *ptr = toupper(*ptr); + ++ptr; + } +} + +/* Return lwalloc'ed hexwkb representation for a GBOX */ +static char * +_box2d_to_hexwkb(const GBOX *bbox, int srid) +{ + POINTARRAY *pa = ptarray_construct(0, 0, 2); + POINT4D p; + LWLINE *line; + char *hex; + size_t sz; + + p.x = bbox->xmin; + p.y = bbox->ymin; + ptarray_set_point4d(pa, 0, &p); + p.x = bbox->xmax; + p.y = bbox->ymax; + ptarray_set_point4d(pa, 1, &p); + line = lwline_construct(srid, NULL, pa); + hex = lwgeom_to_hexwkb( lwline_as_lwgeom(line), WKT_EXTENDED, &sz); + assert(hex[sz-1] == '\0'); + return hex; +} + /* Backend callbacks */ static const char* @@ -159,6 +194,12 @@ cb_loadTopologyByName(const LWT_BE_DATA* be, const char *name) return NULL; } topo->srid = DatumGetInt32(dat); + if ( topo->srid < 0 ) + { + lwnotice("Topology SRID value %d converted to " + "the officially unknown SRID value %d", SRID_UNKNOWN); + topo->srid = SRID_UNKNOWN; + } topo->precision = 0; /* needed ? */ @@ -2472,6 +2513,7 @@ cb_getNodeWithinBox2D ( const LWT_BE_TOPOLOGY* topo, const GBOX* box, int i; int elems_requested = limit; LWT_ISO_NODE* nodes; + char *hexbox; initStringInfo(sql); @@ -2481,10 +2523,10 @@ cb_getNodeWithinBox2D ( const LWT_BE_TOPOLOGY* topo, const GBOX* box, appendStringInfoString(sql, "SELECT "); addNodeFields(sql, fields); } - appendStringInfo(sql, " FROM \"%s\".node WHERE geom && " - "ST_SetSRID(ST_MakeEnvelope(%g,%g,%g,%g),%d)", - topo->name, box->xmin, box->ymin, - box->xmax, box->ymax, topo->srid); + hexbox = _box2d_to_hexwkb(box, topo->srid); + appendStringInfo(sql, " FROM \"%s\".node WHERE geom && '%s'::geometry", + topo->name, hexbox); + lwfree(hexbox); if ( elems_requested == -1 ) { appendStringInfoString(sql, ")"); } else if ( elems_requested > 0 ) { @@ -2542,6 +2584,7 @@ cb_getEdgeWithinBox2D ( const LWT_BE_TOPOLOGY* topo, const GBOX* box, int i; int elems_requested = limit; LWT_ISO_EDGE* edges; + char *hexbox; initStringInfo(sql); @@ -2551,10 +2594,10 @@ cb_getEdgeWithinBox2D ( const LWT_BE_TOPOLOGY* topo, const GBOX* box, appendStringInfoString(sql, "SELECT "); addEdgeFields(sql, fields, 0); } - appendStringInfo(sql, " FROM \"%s\".edge WHERE geom && " - "ST_SetSRID(ST_MakeEnvelope(%g,%g,%g,%g),%d)", - topo->name, box->xmin, box->ymin, - box->xmax, box->ymax, topo->srid); + hexbox = _box2d_to_hexwkb(box, topo->srid); + appendStringInfo(sql, " FROM \"%s\".edge WHERE geom && '%s'::geometry", + topo->name, hexbox); + lwfree(hexbox); if ( elems_requested == -1 ) { appendStringInfoString(sql, ")"); } else if ( elems_requested > 0 ) { @@ -3248,8 +3291,8 @@ Datum ST_GetFaceEdges(PG_FUNCTION_ARGS) funcctx->user_fctx = state; /* - * Build a tuple description for an - * geometry_dump tuple + * Build a tuple description for a + * getfaceedges_returntype tuple */ tupdesc = RelationNameGetTupleDesc("topology.getfaceedges_returntype"); @@ -3946,14 +3989,7 @@ Datum TopoGeo_AddPoint(PG_FUNCTION_ARGS) pt = lwgeom_as_lwpoint(lwgeom); if ( ! pt ) {{ char buf[32]; - char *ptr; - snprintf(buf, 32, "%s", lwtype_name(lwgeom_get_type(lwgeom))); - buf[31] = '\0'; - ptr = buf; - while (*ptr) { - *ptr = toupper(*ptr); - ++ptr; - } + _lwtype_upper_name(lwgeom_get_type(lwgeom), buf, 32); lwgeom_free(lwgeom); PG_FREE_IF_COPY(geom, 1); lwpgerror("Invalid geometry type (%s) passed to TopoGeo_AddPoint" @@ -4004,3 +4040,128 @@ Datum TopoGeo_AddPoint(PG_FUNCTION_ARGS) SPI_finish(); PG_RETURN_INT32(node_id); } + +/* TopoGeo_AddLine(atopology, point, tolerance) */ +Datum TopoGeo_AddLine(PG_FUNCTION_ARGS); +PG_FUNCTION_INFO_V1(TopoGeo_AddLine); +Datum TopoGeo_AddLine(PG_FUNCTION_ARGS) +{ + text* toponame_text; + char* toponame; + double tol; + LWT_ELEMID *elems; + int nelems; + GSERIALIZED *geom; + LWGEOM *lwgeom; + LWLINE *ln; + LWT_TOPOLOGY *topo; + FuncCallContext *funcctx; + MemoryContext oldcontext, newcontext; + FACEEDGESSTATE *state; + Datum result; + LWT_ELEMID id; + + if (SRF_IS_FIRSTCALL()) + { + POSTGIS_DEBUG(1, "TopoGeo_AddLine first call"); + funcctx = SRF_FIRSTCALL_INIT(); + newcontext = funcctx->multi_call_memory_ctx; + + if ( PG_ARGISNULL(0) || PG_ARGISNULL(1) || PG_ARGISNULL(2) ) { + lwpgerror("SQL/MM Spatial exception - null argument"); + PG_RETURN_NULL(); + } + + toponame_text = PG_GETARG_TEXT_P(0); + toponame = text2cstring(toponame_text); + PG_FREE_IF_COPY(toponame_text, 0); + + geom = PG_GETARG_GSERIALIZED_P(1); + lwgeom = lwgeom_from_gserialized(geom); + ln = lwgeom_as_lwline(lwgeom); + if ( ! ln ) {{ + char buf[32]; + _lwtype_upper_name(lwgeom_get_type(lwgeom), buf, 32); + lwgeom_free(lwgeom); + PG_FREE_IF_COPY(geom, 1); + lwpgerror("Invalid geometry type (%s) passed to " + "TopoGeo_AddLinestring, expected LINESTRING", buf); + PG_RETURN_NULL(); + }} + + tol = PG_GETARG_FLOAT8(2); + if ( tol < 0 ) + { + PG_FREE_IF_COPY(geom, 1); + lwpgerror("Tolerance must be >=0"); + PG_RETURN_NULL(); + } + + if ( SPI_OK_CONNECT != SPI_connect() ) { + lwpgerror("Could not connect to SPI"); + PG_RETURN_NULL(); + } + be_data.data_changed = false; + + { + int pre = be_data.topoLoadFailMessageFlavor; + be_data.topoLoadFailMessageFlavor = 1; + topo = lwt_LoadTopology(be_iface, toponame); + be_data.topoLoadFailMessageFlavor = pre; + } + oldcontext = MemoryContextSwitchTo( newcontext ); + pfree(toponame); + if ( ! topo ) { + /* should never reach this point, as lwerror would raise an exception */ + SPI_finish(); + PG_RETURN_NULL(); + } + + POSTGIS_DEBUG(1, "Calling lwt_AddLine"); + elems = lwt_AddLine(topo, ln, tol, &nelems); + POSTGIS_DEBUG(1, "lwt_AddLine returned"); + lwgeom_free(lwgeom); + PG_FREE_IF_COPY(geom, 1); + lwt_FreeTopology(topo); + + if ( nelems < 0 ) { + /* should never reach this point, as lwerror would raise an exception */ + SPI_finish(); + PG_RETURN_NULL(); + } + + state = lwalloc(sizeof(FACEEDGESSTATE)); + state->elems = elems; + state->nelems = nelems; + state->curr = 0; + funcctx->user_fctx = state; + + POSTGIS_DEBUG(1, "TopoGeo_AddLinestring calling SPI_finish"); + + MemoryContextSwitchTo(oldcontext); + + SPI_finish(); + } + + POSTGIS_DEBUG(1, "Per-call invocation"); + + /* stuff done on every call of the function */ + funcctx = SRF_PERCALL_SETUP(); + + /* get state */ + state = funcctx->user_fctx; + + if ( state->curr == state->nelems ) + { + POSTGIS_DEBUG(1, "We're done, cleaning up all"); + SRF_RETURN_DONE(funcctx); + } + + id = state->elems[state->curr++]; + POSTGIS_DEBUGF(1, "TopoGeo_AddLinestring: cur:%d, val:" INT64_FORMAT, + state->curr-1, id); + + result = Int32GetDatum((int32)id); + + SRF_RETURN_NEXT(funcctx, result); +} diff --git a/topology/sql/populate.sql.in b/topology/sql/populate.sql.in index 0619960f6..2fef39c8a 100644 --- a/topology/sql/populate.sql.in +++ b/topology/sql/populate.sql.in @@ -679,206 +679,8 @@ CREATE OR REPLACE FUNCTION topology.TopoGeo_AddPoint(atopology varchar, apoint g -- }{ CREATE OR REPLACE FUNCTION topology.TopoGeo_addLinestring(atopology varchar, aline geometry, tolerance float8 DEFAULT 0) RETURNS SETOF int AS -$$ -DECLARE - rec RECORD; - rec2 RECORD; - sql TEXT; - set1 GEOMETRY; - set2 GEOMETRY; - snapped GEOMETRY; - noded GEOMETRY; - start_node INTEGER; - end_node INTEGER; - id INTEGER; - inodes GEOMETRY; - iedges GEOMETRY; - tol float8; -BEGIN - - -- 0. Check arguments - IF geometrytype(aline) != 'LINESTRING' THEN - RAISE EXCEPTION 'Invalid geometry type (%) passed to TopoGeo_AddLinestring, expected LINESTRING', geometrytype(aline); - END IF; - - -- Get tolerance, if 0 was given - tol := COALESCE( NULLIF(tolerance, 0), topology._st_mintolerance(atopology, aline) ); - - -- 1. Self-node - noded := ST_UnaryUnion(aline); -#ifdef POSTGIS_TOPOLOGY_DEBUG - RAISE DEBUG 'Self-noded: %', ST_AsText(noded); -#endif - - -- 2. Node to edges falling within tol distance - sql := 'WITH nearby AS ( SELECT e.geom FROM ' - || quote_ident(atopology) - || '.edge e WHERE ST_DWithin(e.geom, $1,' - || tol || ') ) SELECT st_collect(geom) FROM nearby'; -#ifdef POSTGIS_TOPOLOGY_DEBUG - RAISE DEBUG '%', sql; -#endif - EXECUTE sql INTO iedges USING noded; - IF iedges IS NOT NULL THEN - -#ifdef POSTGIS_TOPOLOGY_DEBUG - RAISE DEBUG 'Intersecting edges: %', ST_AsText(iedges); -#endif - - snapped := ST_Snap(noded, iedges, tol); -#ifdef POSTGIS_TOPOLOGY_DEBUG - RAISE DEBUG 'Snapped to edges: %', ST_AsText(snapped); -#endif - - noded := ST_Difference(snapped, iedges); -#ifdef POSTGIS_TOPOLOGY_DEBUG - RAISE DEBUG 'Difference: %', ST_AsText(noded); -#endif - - set1 := ST_Intersection(snapped, iedges); -#ifdef POSTGIS_TOPOLOGY_DEBUG - RAISE DEBUG 'Intersection: %', ST_AsText(set1); -#endif - - set2 := ST_LineMerge(set1); -#ifdef POSTGIS_TOPOLOGY_DEBUG - RAISE DEBUG 'LineMerged intersection: %', ST_AsText(set2); -#endif - - noded := ST_Union(noded, set2); -#ifdef POSTGIS_TOPOLOGY_DEBUG - RAISE DEBUG 'Unioned: %', ST_AsText(noded); -#endif - - END IF; - - -- 2.1. Node with existing nodes within tol - -- TODO: check if we should be only considering _isolated_ nodes! - sql := 'WITH nearby AS ( SELECT n.geom FROM ' - || quote_ident(atopology) - || '.node n WHERE ST_DWithin(n.geom, $1, ' - || tol || ') ) SELECT st_collect(geom) FROM nearby;'; -#ifdef POSTGIS_TOPOLOGY_DEBUG - RAISE DEBUG '%', sql; -#endif - EXECUTE sql INTO inodes USING noded; - - IF inodes IS NOT NULL THEN -- { -#ifdef POSTGIS_TOPOLOGY_DEBUG - RAISE DEBUG 'Intersecting nodes: %', ST_AsText(inodes); -#endif - - -- TODO: consider snapping once against all elements - --- (rather than once with edges and once with nodes) - noded := ST_Snap(noded, inodes, tol); -#ifdef POSTGIS_TOPOLOGY_DEBUG - RAISE DEBUG 'Snapped to nodes: %', ST_AsText(noded); -#endif - - FOR rec IN SELECT (ST_Dump(inodes)).geom - LOOP - -- Use the node to split edges - SELECT ST_Collect(geom) - FROM ST_Dump(ST_Split(noded, rec.geom)) - INTO STRICT noded; -#ifdef POSTGIS_TOPOLOGY_DEBUG - RAISE DEBUG 'Split by %: %', ST_AsText(rec.geom), ST_AsText(noded); -#endif - END LOOP; -#ifdef POSTGIS_TOPOLOGY_DEBUG - RAISE DEBUG 'Split: %', ST_AsText(noded); -#endif - - -- 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 - noded := ST_UnaryUnion(noded); -#ifdef POSTGIS_TOPOLOGY_DEBUG - RAISE DEBUG 'Self-unioned again: %', ST_AsText(noded); -#endif - END IF; -- } - - -- 3. For each (now-noded) segment, insert an edge - FOR rec IN SELECT (ST_Dump(noded)).geom LOOP - - -- TODO: skip point elements ? - -#ifdef POSTGIS_TOPOLOGY_DEBUG - RAISE DEBUG 'Adding edge %', ST_AsText(rec.geom); -#endif - - start_node := topology.TopoGeo_AddPoint(atopology, - ST_StartPoint(rec.geom), - tol); -#ifdef POSTGIS_TOPOLOGY_DEBUG - RAISE DEBUG ' Start Node: %', start_node; -#endif - - end_node := topology.TopoGeo_AddPoint(atopology, - ST_EndPoint(rec.geom), - tol); -#ifdef POSTGIS_TOPOLOGY_DEBUG - RAISE DEBUG ' End Node: %', end_node; -#endif - - -- Added endpoints may have drifted due to tolerance, so - -- we need to re-snap the edge to the new nodes before adding it - sql := 'SELECT n1.geom as sn, n2.geom as en FROM ' || quote_ident(atopology) - || '.node n1, ' || quote_ident(atopology) - || '.node n2 WHERE n1.node_id = ' - || start_node || ' AND n2.node_id = ' || end_node; -#ifdef POSTGIS_TOPOLOGY_DEBUG - RAISE DEBUG '%', sql; -#endif - - EXECUTE sql INTO STRICT rec2; - - snapped := ST_SetPoint( - ST_SetPoint(rec.geom, ST_NPoints(rec.geom)-1, rec2.en), - 0, rec2.sn); - - /* We might have introduced an invalidity (TODO: check this out) */ - snapped := ST_CollectionExtract(ST_MakeValid(snapped), 2); -#ifdef POSTGIS_TOPOLOGY_DEBUG - RAISE DEBUG 'Cleaned edge: %', ST_AsText(snapped); -#endif - - - -- Check if the so-snapped edge collapsed (see #1650) - IF ST_IsEmpty(snapped) THEN -#ifdef POSTGIS_TOPOLOGY_DEBUG - RAISE DEBUG 'Edge collapsed'; -#endif - CONTINUE; - END IF; - - -- Check if the so-snapped edge _now_ exists - sql := 'SELECT edge_id FROM ' || quote_ident(atopology) - || '.edge_data WHERE ST_Equals(geom, $1)'; -#ifdef POSTGIS_TOPOLOGY_DEBUG - RAISE DEBUG '%', sql; -#endif - EXECUTE sql INTO id USING snapped; - IF id IS NULL THEN - id := topology.ST_AddEdgeModFace(atopology, start_node, end_node, - snapped); -#ifdef POSTGIS_TOPOLOGY_DEBUG - RAISE DEBUG 'New edge id: %', id; -#endif - ELSE -#ifdef POSTGIS_TOPOLOGY_DEBUG - RAISE DEBUG 'Old edge id: %', id; -#endif - END IF; - - RETURN NEXT id; - - END LOOP; - - RETURN; -END -$$ -LANGUAGE 'plpgsql'; + 'MODULE_PATHNAME', 'TopoGeo_AddLine' + LANGUAGE 'c' VOLATILE; --} TopoGeo_addLinestring --{