[postgis-tickets] [SCM] PostGIS branch main updated. 3.1.0rc1-394-gb565cc6
git at osgeo.org
git at osgeo.org
Wed Aug 11 07:28:29 PDT 2021
This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "PostGIS".
The branch, main has been updated
via b565cc631df71e8e8081a747e5dd211c486e2f35 (commit)
via 6f9afe7e25c561ebcc0cdc009d305133cf2bb279 (commit)
via ec3f582abbf9d22d07338521317f8ab94d144c29 (commit)
from d3f5c2e14c450ee0ca8f7f8dbd30f1f1b8ce1ae0 (commit)
Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.
- Log -----------------------------------------------------------------
commit b565cc631df71e8e8081a747e5dd211c486e2f35
Author: Sandro Santilli <strk at kbt.io>
Date: Tue Aug 10 13:58:18 2021 +0200
Rewrite GetFaceContainingPoint in C
Closes #4962
diff --git a/liblwgeom/liblwgeom_topo.h b/liblwgeom/liblwgeom_topo.h
index 665e9c5..c49ba0b 100644
--- a/liblwgeom/liblwgeom_topo.h
+++ b/liblwgeom/liblwgeom_topo.h
@@ -18,7 +18,7 @@
*
**********************************************************************
*
- * Copyright (C) 2015-2017 Sandro Santilli <strk at kbt.io>
+ * Copyright (C) 2015-2021 Sandro Santilli <strk at kbt.io>
*
**********************************************************************/
@@ -357,20 +357,6 @@ typedef struct LWT_BE_CALLBACKS_T {
LWT_ISO_FACE *(*getFaceById)(const LWT_BE_TOPOLOGY *topo, const LWT_ELEMID *ids, uint64_t *numelems, int fields);
/**
- * Get face containing point
- *
- * @param topo the topology to act upon
- * @param pt the query point
- *
- * @return a face identifier (0 if no face contains the point)
- * or -1 on error (@see lastErrorMessage)
- */
- LWT_ELEMID (*getFaceContainingPoint) (
- const LWT_BE_TOPOLOGY* topo,
- const LWPOINT* pt
- );
-
- /**
* Update TopoGeometry objects after an edge split event
*
* @param topo the topology to act upon
@@ -846,6 +832,27 @@ typedef struct LWT_BE_CALLBACKS_T {
LWT_ELEMID rem_edge
);
+
+ /**
+ * Get closest edge to a given point
+ *
+ * @param topo the topology to act upon
+ * @param pt the query point
+ * @param numelems output parameter, gets number of elements found
+ * or UINT64_MAX on error (@see lastErrorMessage)
+ * @param fields fields to be filled in the returned structure, see
+ * LWT_COL_EDGE_* macros
+ *
+ * @return an array of 1 edges or null in the following cases:
+ * - no edges are in the topology ("numelems" is set to 0)
+ * - error ("numelems" is set to UINT64_MAX)
+ *
+ */
+ LWT_ISO_EDGE *(*getClosestEdge)(const LWT_BE_TOPOLOGY *topo,
+ const LWPOINT *pt,
+ uint64_t *numelems,
+ int fields);
+
} LWT_BE_CALLBACKS;
@@ -993,7 +1000,20 @@ LWT_ELEMID lwt_GetEdgeByPoint(LWT_TOPOLOGY *topo, LWPOINT *pt, double tol);
* or edge).
* The liblwgeom error handler will be invoked in case of error.
*/
-LWT_ELEMID lwt_GetFaceByPoint(LWT_TOPOLOGY *topo, LWPOINT *pt, double tol);
+LWT_ELEMID lwt_GetFaceByPoint(LWT_TOPOLOGY *topo, const LWPOINT *pt, double tol);
+
+/**
+ * Find the face-id of the face properly containing a given point
+ *
+ * @param topo the topology to operate on
+ * @param point the point to use for query
+ *
+ * @return a face identifier if one is found (0 if universe), -1
+ * on error (point intersects non-dangling edge).
+ * The liblwgeom error handler will be invoked in case of error.
+ */
+LWT_ELEMID lwt_GetFaceContainingPoint(LWT_TOPOLOGY* topo, const LWPOINT* pt);
+
/*******************************************************************
diff --git a/liblwgeom/liblwgeom_topo_internal.h b/liblwgeom/liblwgeom_topo_internal.h
index de26433..cbd7ade 100644
--- a/liblwgeom/liblwgeom_topo_internal.h
+++ b/liblwgeom/liblwgeom_topo_internal.h
@@ -66,7 +66,7 @@ int lwt_be_ExistsEdgeIntersectingPoint(LWT_TOPOLOGY* topo, LWPOINT* pt);
LWT_ELEMID lwt_be_getNextEdgeId(LWT_TOPOLOGY* topo);
LWT_ISO_EDGE *lwt_be_getEdgeById(LWT_TOPOLOGY *topo, const LWT_ELEMID *ids, uint64_t *numelems, int fields);
LWT_ISO_EDGE *lwt_be_getEdgeWithinDistance2D(LWT_TOPOLOGY *topo,
- LWPOINT *pt,
+ const LWPOINT *pt,
double dist,
uint64_t *numelems,
int fields,
@@ -77,8 +77,6 @@ lwt_be_updateEdges(LWT_TOPOLOGY* topo, const LWT_ISO_EDGE* sel_edge, int sel_fie
int
lwt_be_deleteEdges(LWT_TOPOLOGY* topo, const LWT_ISO_EDGE* sel_edge, int sel_fields);
-LWT_ELEMID lwt_be_getFaceContainingPoint(LWT_TOPOLOGY* topo, LWPOINT* pt);
-
int lwt_be_updateTopoGeomEdgeSplit(LWT_TOPOLOGY* topo, LWT_ELEMID split_edge, LWT_ELEMID new_edge1, LWT_ELEMID new_edge2);
diff --git a/liblwgeom/lwgeom_topo.c b/liblwgeom/lwgeom_topo.c
index 8b7e2cd..427aad8 100644
--- a/liblwgeom/lwgeom_topo.c
+++ b/liblwgeom/lwgeom_topo.c
@@ -248,7 +248,7 @@ lwt_be_getNodeByFace(LWT_TOPOLOGY *topo, const LWT_ELEMID *ids, uint64_t *numele
LWT_ISO_EDGE *
lwt_be_getEdgeWithinDistance2D(LWT_TOPOLOGY *topo,
- LWPOINT *pt,
+ const LWPOINT *pt,
double dist,
uint64_t *numelems,
int fields,
@@ -319,13 +319,6 @@ lwt_be_deleteEdges(LWT_TOPOLOGY* topo,
CBT2(topo, deleteEdges, sel_edge, sel_fields);
}
-LWT_ELEMID
-lwt_be_getFaceContainingPoint(LWT_TOPOLOGY* topo, LWPOINT* pt)
-{
- CBT1(topo, getFaceContainingPoint, pt);
-}
-
-
int
lwt_be_updateTopoGeomEdgeSplit(LWT_TOPOLOGY* topo, LWT_ELEMID split_edge, LWT_ELEMID new_edge1, LWT_ELEMID new_edge2)
{
@@ -387,9 +380,6 @@ lwt_be_getRingEdges(LWT_TOPOLOGY *topo, LWT_ELEMID edge, uint64_t *numedges, uin
CBT3(topo, getRingEdges, edge, numedges, limit);
}
-
-/* wrappers of backend wrappers... */
-
int
lwt_be_ExistsCoincidentNode(LWT_TOPOLOGY* topo, LWPOINT* pt)
{
@@ -416,6 +406,12 @@ lwt_be_ExistsEdgeIntersectingPoint(LWT_TOPOLOGY* topo, LWPOINT* pt)
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);
+}
+
/************************************************************************
*
* Utility functions
@@ -554,7 +550,7 @@ _lwt_AddIsoNode( LWT_TOPOLOGY* topo, LWT_ELEMID face,
if ( checkFace && ( face == -1 || ! skipISOChecks ) )
{
- foundInFace = lwt_be_getFaceContainingPoint(topo, pt); /*x*/
+ foundInFace = lwt_GetFaceContainingPoint(topo, pt); /*x*/
if ( foundInFace == -1 ) {
lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
return -1;
@@ -3632,7 +3628,7 @@ lwt_MoveIsoNode(LWT_TOPOLOGY* topo, LWT_ELEMID nid, LWPOINT *pt)
/* Check that the new point is in the same containing face !
* See https://trac.osgeo.org/postgis/ticket/3232 */
- newPointFace = lwt_be_getFaceContainingPoint(topo, pt);
+ newPointFace = lwt_GetFaceContainingPoint(topo, pt);
if ( newPointFace == -1 ) {
lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
return -1;
@@ -4803,7 +4799,7 @@ lwt_GetEdgeByPoint(LWT_TOPOLOGY *topo, LWPOINT *pt, double tol)
}
LWT_ELEMID
-lwt_GetFaceByPoint(LWT_TOPOLOGY *topo, LWPOINT *pt, double tol)
+lwt_GetFaceByPoint(LWT_TOPOLOGY *topo, const LWPOINT *pt, double tol)
{
LWT_ELEMID id = 0;
LWT_ISO_EDGE *elem;
@@ -4814,7 +4810,7 @@ lwt_GetFaceByPoint(LWT_TOPOLOGY *topo, LWPOINT *pt, double tol)
LWT_COL_EDGE_FACE_RIGHT;
LWGEOM *qp = lwpoint_as_lwgeom(pt);
- id = lwt_be_getFaceContainingPoint(topo, pt);
+ id = lwt_GetFaceContainingPoint(topo, pt);
if ( id == -1 ) {
lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_iface));
return -1;
@@ -7054,3 +7050,238 @@ lwt_Polygonize(LWT_TOPOLOGY* topo)
return 0;
}
+
+LWT_ELEMID
+lwt_GetFaceContainingPoint(LWT_TOPOLOGY* topo, const LWPOINT* pt)
+{
+ LWT_ISO_EDGE* closestEdge;
+ LWT_ISO_EDGE* edges;
+ uint64_t numedges;
+ LWGEOM *shortestLine;
+ const POINT2D *shortestLineP0, *shortestLineP1;
+ int closestSegmentIndex;
+ int closestSegmentSide;
+ const POINT2D *closestSegmentP0, *closestSegmentP1;
+ double shortestLineAzimuth;
+ LWPOINT *closestPointOnEdge;
+ LWPOINT *closestEdgeEndpoint;
+ LWT_ELEMID closestNode = 0;
+ double dist;
+ int containingFace = -1;
+ int i;
+
+ closestEdge = lwt_be_getClosestEdge( topo, pt, &numedges,
+ LWT_COL_EDGE_GEOM|
+ LWT_COL_EDGE_EDGE_ID|
+ LWT_COL_EDGE_START_NODE|
+ LWT_COL_EDGE_END_NODE|
+ LWT_COL_EDGE_FACE_LEFT|
+ LWT_COL_EDGE_FACE_RIGHT
+ );
+ if (numedges == UINT64_MAX)
+ {
+ lwerror("Backend error: %s", lwt_be_lastErrorMessage(topo->be_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;
+ }
+
+
+ /* Compute shortest line from query point to closest edge */
+
+ shortestLine = lwgeom_closest_line(lwline_as_lwgeom(closestEdge->geom), lwpoint_as_lwgeom(pt));
+
+ LWDEBUGGF(1, shortestLine, "Shortest line to closest edge %d", closestEdge->edge_id);
+
+ /* Check if the closest point on the edge is an edge endpoint */
+
+ closestPointOnEdge = lwline_get_lwpoint(lwgeom_as_lwline(shortestLine), 0);
+
+ LWDEBUGG(1, lwpoint_as_lwgeom(closestPointOnEdge), "Closest point on closest edge");
+
+ closestEdgeEndpoint = lwline_get_lwpoint(closestEdge->geom, 0);
+ if ( lwpoint_same(closestEdgeEndpoint, closestPointOnEdge) ) {
+ closestNode = closestEdge->start_node;
+ } else if ( ! lwline_is_closed(closestEdge->geom) ) {
+ lwpoint_free(closestEdgeEndpoint);
+ closestEdgeEndpoint = lwline_get_lwpoint(closestEdge->geom, closestEdge->geom->points->npoints - 1);
+ if ( lwpoint_same(closestEdgeEndpoint, closestPointOnEdge) ) {
+ closestNode = closestEdge->end_node;
+ }
+ }
+ lwpoint_free(closestEdgeEndpoint);
+ lwpoint_free(closestPointOnEdge);
+
+ if ( closestNode != 0 )
+ {
+ LWDEBUGF(1, "Closest point is node %d", closestNode);
+ dist = lwgeom_length(shortestLine);
+ if ( dist == 0 )
+ {
+ lwgeom_free(shortestLine);
+ 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; i<numedges; ++i)
+ {
+ if ( edges[i].face_left != containingFace ||
+ edges[i].face_right != containingFace )
+ {
+ _lwt_release_edges(edges, numedges);
+ _lwt_release_edges(closestEdge, 1);
+ lwerror("Two or more faces found");
+ return -1;
+ }
+ }
+ if (numedges < 1 )
+ {
+ lwerror("Unexpected backend return: getEdgeByNode(%d) returns no edges when we previously found edge %d ending on that node",
+ closestNode, closestEdge->edge_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 */
+ shortestLineP0 = getPoint2d_cp(((LWLINE *)shortestLine)->points, 0);
+ shortestLineP1 = getPoint2d_cp(((LWLINE *)shortestLine)->points, 1);
+ if ( ! azimuth_pt_pt(shortestLineP0, shortestLineP1, &shortestLineAzimuth) ) {
+ lwerror("error computing azimuth of shortestLine [%.15g %.15g,%.15g %.15g]",
+ shortestLineP0->x, shortestLineP0->y,
+ shortestLineP1->x, shortestLineP1->y);
+ return -1;
+ }
+
+ LWDEBUGF(1, "ShortestLine azimuth is %g", shortestLineAzimuth);
+
+ edgeend ee;
+ ee.myaz = shortestLineAzimuth;
+ 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;
+
+ }
+
+ LWDEBUGF(1, "Closest point is NOT a node", closestNode);
+
+ /* 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;
+ }
+
+ dist = lwgeom_length(shortestLine);
+ if ( dist == 0 )
+ {
+ /* We checked the dangling case above */
+ _lwt_release_edges(closestEdge, 1);
+ lwerror("Two or more faces found");
+ return -1;
+ }
+
+ /* find closest segment on closestEdge and use lw_segment_side
+ * to determine on which side our query point falls */
+
+ shortestLineP0 = getPoint2d_cp(((LWLINE *)shortestLine)->points, 1);
+ closestSegmentIndex = ptarray_closest_segment_2d(closestEdge->geom->points, shortestLineP0, NULL);
+ LWDEBUGF(1, "Closest segment to edge %d is %d", closestEdge->edge_id, closestSegmentIndex);
+ closestSegmentP0 = getPoint2d_cp(closestEdge->geom->points, closestSegmentIndex);
+ closestSegmentP1 = getPoint2d_cp(closestEdge->geom->points, ++closestSegmentIndex);
+ LWDEBUGF(1, "Closest segment to edge %d is LINESTRING(%g %g, %g %g)",
+ closestEdge->edge_id,
+ closestSegmentP0->x,
+ closestSegmentP0->y,
+ closestSegmentP1->x,
+ closestSegmentP1->y
+ );
+
+ /* Find on which side of the segment the query point lays */
+
+ do {
+ containingFace = -1;
+
+ shortestLineP1 = getPoint2d_cp(((LWLINE *)shortestLine)->points, 1);
+ LWDEBUGF(1, "Calling lw_segment_side((%g,%g),(%g,%g),(%g,%g)",
+ closestSegmentP0->x,
+ closestSegmentP0->y,
+ closestSegmentP1->x,
+ closestSegmentP1->y,
+ shortestLineP1->x,
+ shortestLineP1->y
+ );
+
+ closestSegmentSide = lw_segment_side(closestSegmentP0, closestSegmentP1, shortestLineP1);
+ LWDEBUGF(1, "Side of closest segment query point falls on: %d", closestSegmentSide);
+
+ if ( closestSegmentSide == -1 ) /* left */
+ {
+ containingFace = closestEdge->face_left;
+ break;
+ }
+ else if ( closestSegmentSide == 1 ) /* right */
+ {
+ containingFace = closestEdge->face_right;
+ break;
+ }
+
+ /* If the point was collinear to the closest segment
+ * we need to check the next segment, as if this was
+ * the last one we would have entered a previous branch
+ * selecting the case in which the closest point is
+ * a NODE, and ptarray_closest_segment_2d should be
+ * always returning the *first* segment so there must
+ * be a next one.
+ */
+ closestSegmentP0 = closestSegmentP1;
+ closestSegmentP1 = getPoint2d_cp(closestEdge->geom->points, ++closestSegmentIndex);
+
+ } while (1);
+
+ _lwt_release_edges(closestEdge, 1);
+ lwgeom_free(shortestLine);
+ return containingFace;
+}
diff --git a/topology/postgis_topology.c b/topology/postgis_topology.c
index d23831d..9f0756d 100644
--- a/topology/postgis_topology.c
+++ b/topology/postgis_topology.c
@@ -3,7 +3,7 @@
* PostGIS - Spatial Types for PostgreSQL
* http://postgis.net
*
- * Copyright (C) 2015 Sandro Santilli <strk at kbt.io>
+ * Copyright (C) 2015-2021 Sandro Santilli <strk at kbt.io>
*
* This is free software; you can redistribute and/or modify it under
* the terms of the GNU General Public Licence. See the COPYING file.
@@ -2976,66 +2976,85 @@ cb_updateTopoGeomEdgeHeal ( const LWT_BE_TOPOLOGY* topo,
return 1;
}
-static LWT_ELEMID
-cb_getFaceContainingPoint( const LWT_BE_TOPOLOGY* topo, const LWPOINT* pt )
+/* Return pointer to closest edge.
+ * NULL if no edges are in the topology (*numelems=0).
+ * or there's an error (*numedges=UINT64_MAX)
+ */
+static LWT_ISO_EDGE *
+cb_getClosestEdge( const LWT_BE_TOPOLOGY* topo, const LWPOINT* pt, uint64_t *numedges, int fields )
{
MemoryContext oldcontext = CurrentMemoryContext;
+ HeapTuple row;
int spi_result;
StringInfoData sqldata;
StringInfo sql = &sqldata;
- bool isnull;
- Datum dat;
- LWT_ELEMID face_id;
GSERIALIZED *pts;
- Datum values[2];
- Oid argtypes[2];
- text *toponameText;
-
- initStringInfo(sql);
+ Datum values[1];
+ Oid argtypes[1];
+ LWT_ISO_EDGE* edges;
pts = geometry_serialize(lwpoint_as_lwgeom(pt));
if ( ! pts )
{
cberror(topo->be_data, "%s:%d: could not serialize query point",
__FILE__, __LINE__);
- return -1;
+ *numedges = UINT64_MAX;
+ return NULL;
}
- appendStringInfo(sql, "SELECT topology.GetFaceContainingPoint($1, $2)");
- toponameText = cstring_to_text(topo->name);
- values[0] = PointerGetDatum(toponameText);
- argtypes[0] = TEXTOID;
+ initStringInfo(sql);
+
- values[1] = PointerGetDatum(pts);
- argtypes[1] = topo->geometryOID;
+ /* TODO: avoid executing ST_ShortestLine twice (to be confirmed it
+ * if it is called twice
+ */
+ appendStringInfoString(sql, "SELECT ");
+ addEdgeFields(sql, fields, 0);
+ appendStringInfo(sql, " FROM \"%s\".edge_data ORDER BY ", topo->name);
+
+#if POSTGIS_PGSQL_VERSION < 95
+ appendStringInfo(sql, "ST_Length(ST_ShortestLine(geom, $1))");
+#else
+ appendStringInfo(sql, "geom <-> $1");
+#endif
- spi_result = SPI_execute_with_args(sql->data, 2, argtypes, values, NULL,
+ appendStringInfo(sql, " ASC, edge_id ASC LIMIT 1");
+
+ POSTGIS_DEBUGF(1, "cb_getClosestEdge query: %s", sql->data);
+
+ values[0] = PointerGetDatum(pts);
+ argtypes[0] = topo->geometryOID;
+
+ spi_result = SPI_execute_with_args(sql->data, 1, argtypes, values, NULL,
!topo->be_data->data_changed, 1);
MemoryContextSwitchTo( oldcontext ); /* switch back */
pfree(pts); /* not needed anymore */
- if ( spi_result != SPI_OK_SELECT || SPI_processed != 1 )
+ if ( spi_result != SPI_OK_SELECT )
{
cberror(topo->be_data, "unexpected return (%d) from query execution: %s",
spi_result, sql->data);
pfree(sqldata.data);
- return -1;
+ *numedges = UINT64_MAX;
+ return NULL;
}
- dat = SPI_getbinval( SPI_tuptable->vals[0],
- SPI_tuptable->tupdesc, 1, &isnull );
- if ( isnull )
+ if ( SPI_processed == 0 )
{
- SPI_freetuptable(SPI_tuptable);
- cberror(topo->be_data, "unexpected return (%d) from query execution: %s",
- spi_result, sql->data);
+ /* No edges in topology, point is in universal face */
pfree(sqldata.data);
- return -1;
+ *numedges = 0;
+ return NULL;
}
- pfree(sqldata.data);
- face_id = DatumGetInt32(dat);
+ /* Got closest edge, return it */
+ *numedges = 1;
+
+ edges = palloc( sizeof(LWT_ISO_EDGE) );
+ row = SPI_tuptable->vals[0];
+ fillEdgeFields(edges, row, SPI_tuptable->tupdesc, fields);
SPI_freetuptable(SPI_tuptable);
- return face_id;
+
+ return edges;
}
static int
@@ -3389,7 +3408,6 @@ static LWT_BE_CALLBACKS be_callbacks =
cb_insertEdges,
cb_updateEdges,
cb_getFacesById,
- cb_getFaceContainingPoint,
cb_updateTopoGeomEdgeSplit,
cb_deleteEdges,
cb_getNodeWithinBox2D,
@@ -3415,7 +3433,8 @@ static LWT_BE_CALLBACKS be_callbacks =
cb_updateTopoGeomEdgeHeal,
cb_getFaceWithinBox2D,
cb_checkTopoGeomRemIsoNode,
- cb_checkTopoGeomRemIsoEdge
+ cb_checkTopoGeomRemIsoEdge,
+ cb_getClosestEdge
};
static void
@@ -5314,3 +5333,64 @@ Datum GetRingEdges(PG_FUNCTION_ARGS)
SRF_RETURN_NEXT(funcctx, result);
}
+
+/* GetFaceContainingPoint(atopology, point) */
+Datum GetFaceContainingPoint(PG_FUNCTION_ARGS);
+PG_FUNCTION_INFO_V1(GetFaceContainingPoint);
+Datum GetFaceContainingPoint(PG_FUNCTION_ARGS)
+{
+ text* toponame_text;
+ char* toponame;
+ LWT_ELEMID face_id;
+ GSERIALIZED *geom;
+ LWGEOM *lwgeom;
+ LWPOINT *pt;
+ LWT_TOPOLOGY *topo;
+
+ toponame_text = PG_GETARG_TEXT_P(0);
+ toponame = text_to_cstring(toponame_text);
+ PG_FREE_IF_COPY(toponame_text, 0);
+
+ geom = PG_GETARG_GSERIALIZED_P(1);
+ lwgeom = lwgeom_from_gserialized(geom);
+ pt = lwgeom_as_lwpoint(lwgeom);
+ if ( ! pt )
+ {
+ lwgeom_free(lwgeom);
+ PG_FREE_IF_COPY(geom, 1);
+ lwpgerror("Second argument must be a point geometry");
+ PG_RETURN_NULL();
+ }
+
+ if ( SPI_OK_CONNECT != SPI_connect() )
+ {
+ lwpgerror("Could not connect to SPI");
+ PG_RETURN_NULL();
+ }
+
+ topo = lwt_LoadTopology(be_iface, toponame);
+ 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_GetFaceContainingPoint");
+ face_id = lwt_GetFaceContainingPoint(topo, pt);
+ POSTGIS_DEBUG(1, "lwt_GetFaceContainingPoint returned");
+ lwgeom_free(lwgeom);
+ PG_FREE_IF_COPY(geom, 1);
+ lwt_FreeTopology(topo);
+
+ if ( face_id == -1 )
+ {
+ /* should never reach this point, as lwerror would raise an exception */
+ SPI_finish();
+ PG_RETURN_NULL();
+ }
+
+ SPI_finish();
+ PG_RETURN_INT32(face_id);
+}
diff --git a/topology/sql/query/GetFaceContainingPoint.sql.in b/topology/sql/query/GetFaceContainingPoint.sql.in
index 86c6337..bd48bd2 100644
--- a/topology/sql/query/GetFaceContainingPoint.sql.in
+++ b/topology/sql/query/GetFaceContainingPoint.sql.in
@@ -25,287 +25,6 @@ CREATE OR REPLACE FUNCTION topology.GetFaceContainingPoint(
atopology text,
apoint geometry
)
-RETURNS INT
-AS $BODY$
-DECLARE
- closestEdge RECORD;
- closestCWEdgeEnd RECORD;
- rec RECORD;
- closestNode INT;
- azimuth FLOAT8;
- closestVertex GEOMETRY;
- lastDistance FLOAT8;
- ringLeft GEOMETRY;
- ringRight GEOMETRY;
- sql TEXT;
-BEGIN
-
- sql := format(
- $$
- SELECT
- edge_id,
- geom, -- for Start/EndPoint query
- start_node,
- end_node,
- left_face,
- right_face,
- ST_ShortestLine(geom, $1) AS shortestLine,
- ST_Length(ST_ShortestLine(geom, $1)) AS dist
- FROM %1$I.edge
- ORDER BY
-#if POSTGIS_PGSQL_VERSION < 95
- ST_Length(ST_ShortestLine(geom, $1))
-#else
- geom <-> $1
-#endif
- ASC,
- edge_id ASC
- LIMIT 1
- $$,
- atopology
- );
-#ifdef POSTGIS_TOPOLOGY_DEBUG
- RAISE DEBUG 'Query: %', sql;
- RAISE DEBUG 'Point: %', ST_AsEWKT(apoint);
-#endif
- EXECUTE sql
- USING apoint
- INTO closestEdge;
-
- IF closestEdge IS NULL
- THEN
- -- If no edge is found the point
- -- is in the universal face
- RETURN 0;
- END IF;
-
- -- Check if the closest vertex is a node
- closestVertex := ST_StartPoint(closestEdge.shortestLine);
- IF ST_Equals(closestVertex, ST_StartPoint(closestEdge.geom))
- THEN
- closestNode := closestEdge.start_node;
- ELSIF ST_Equals(closestVertex, ST_EndPoint(closestEdge.geom))
- THEN
- closestNode := closestEdge.end_node;
- END IF;
-
-#ifdef POSTGIS_TOPOLOGY_DEBUG
- RAISE DEBUG 'Closest edge is %, shortest line is % (length %)',
- closestEdge.edge_id,
- ST_AsEWKT(closestEdge.shortestLine),
- closestEdge.dist
- ;
-#endif
-
- IF closestNode IS NOT NULL
- THEN --{ Closest point is a node
-
- RAISE DEBUG 'Closest point is node %', closestNode;
-
- IF closestEdge.dist = 0
- THEN --{ 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
-
- -- early exit if the closest edge has
- -- left_face != right_face, don't bother
- -- fetching more edges
- IF closestEdge.left_face != closestEdge.right_face THEN
- RAISE EXCEPTION 'Two or more faces found';
- END IF;
-
- sql := format(
- $$
- SELECT
- left_face,
- right_face
- FROM %1$I.edge_data
- WHERE (
- start_node = $1
- OR end_node = $1
- )
- AND (
- left_face != $2
- OR right_face != $2
- )
- LIMIT 1
- $$,
- atopology
- );
-#ifdef POSTGIS_TOPOLOGY_DEBUG
- RAISE DEBUG 'SQL: %', sql;
-#endif
- EXECUTE sql USING closestNode, closestEdge.left_face INTO rec;
- IF FOUND
- THEN
- RAISE EXCEPTION 'Two or more faces found';
- END IF;
- RETURN closestEdge.left_face;
-
- END IF; --} query point is the node
-
- -- Compute azimuth of segment going from node to query point
- azimuth := ST_Azimuth(
- ST_StartPoint( closestEdge.shortestLine),
- ST_EndPoint(closestEdge.shortestLine)
- );
-
- RAISE DEBUG 'Closest point is disjoint from node %, shortestLine has azimuth %',
- closestNode,
- azimuth;
-
- sql := format(
- $$
- WITH incident_edges AS (
- SELECT
- edge_id,
- start_node,
- end_node,
- left_face,
- right_face,
- ST_RemoveRepeatedPoints(geom) as geom
- FROM %1$I.edge_data
- WHERE start_node = $1
- OR end_node = $1
- ), incident_edges_with_azdiff AS (
- SELECT
- edge_id,
- left_face,
- right_face,
- ST_Azimuth(ST_StartPoint(geom), ST_PointN(geom, 2)) - %2$L as azdiff
- FROM
- incident_edges WHERE start_node = $1
- UNION ALL
- SELECT
- -edge_id,
- left_face,
- right_face,
- ST_Azimuth(ST_EndPoint(geom), ST_PointN(geom, ST_NumPoints(geom)-1)) - %2$L
- FROM incident_edges WHERE end_node = $1
- )
- SELECT
- edge_id,
- left_face,
- right_face,
- CASE WHEN azdiff < 0 THEN
- azdiff + 2 * PI()
- ELSE
- azdiff
- END az
- FROM incident_edges_with_azdiff
- ORDER BY az ASC
- LIMIT 1
- $$,
- atopology,
- azimuth
- );
- RAISE DEBUG 'SQL: %', sql;
- EXECUTE sql USING closestNode INTO closestCWEdgeEnd;
-
- RAISE DEBUG 'Adjiacent edge in CW direction is %, with azimuth %',
- closestCWEdgeEnd.edge_id,
- closestCWEdgeEnd.az
- ;
-
- -- Return the face on the appropriate side of the
- -- most adjiacent edge end
- IF closestCWEdgeEnd.edge_id > 0 THEN --{{
- -- Edge starts on given node, so we pick its left face
- RETURN closestCWEdgeEnd.left_face;
- ELSE -- }{
- -- Edge ends on given node, so we pick its right face
- RETURN closestCWEdgeEnd.right_face;
- END IF; -- }}
-
- END IF; --} closest point is 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.left_face = closestEdge.right_face
- THEN
- RETURN closestEdge.left_face;
- END IF;
-
- -- The query point intersects an edge
- IF closestEdge.dist = 0 THEN
- -- We checked the dangling case above
- RAISE EXCEPTION 'Two or more faces found';
- END IF;
-
- -- TODO: use ST_LineCrossingDirection instead of building rings ?
- -- Beware of https://trac.osgeo.org/postgis/ticket/4935
-
- -- Compute rings on the left and on the right of the closest edge
- sql := format(
- $$
- WITH ring AS (
- SELECT * FROM topology.GetRingEdges(%1$L, $1)
- ), edges AS (
- SELECT
- r.sequence,
- CASE WHEN r.edge > 0 THEN
- e.geom
- ELSE
- ST_Reverse(e.geom)
- END as geom
- FROM
- ring r,
- %1$I.edge e
- WHERE e.edge_id = abs(r.edge)
- )
- SELECT ST_MakePolygon(
- ST_MakeLine(
- geom
- ORDER BY sequence
- )
- )
- FROM edges
- $$,
- atopology
- );
-
- BEGIN
- EXECUTE sql USING closestEdge.edge_id INTO ringLeft;
- EXCEPTION WHEN OTHERS THEN
- IF SQLERRM like '%shell must have at least%' THEN
- RAISE EXCEPTION 'Corrupted topology: ring on the left side of edge % is invalid: %',
- closestEdge.edge_id, SQLERRM;
- ELSE
- RAISE EXCEPTION 'Could not build ring on the left side of edge %: %',
- closestEdge.edge_id, SQLERRM;
- END IF;
- END;
-
- BEGIN
- EXECUTE sql USING -closestEdge.edge_id INTO ringRight;
- EXCEPTION WHEN OTHERS THEN
- IF SQLERRM like '%shell must have at least%' THEN
- RAISE EXCEPTION 'Corrupted topology: ring on the right side of edge % is invalid: %',
- closestEdge.edge_id, SQLERRM;
- ELSE
- RAISE EXCEPTION 'Could not build ring on the right side of edge %: %',
- closestEdge.edge_id, SQLERRM;
- END IF;
- END;
-
- IF ST_IsPolygonCCW(ringLeft) THEN -- ring on left is a shell
- IF ST_Contains(ringLeft, apoint) THEN
- RETURN closestEdge.left_face;
- ELSE
- RETURN closestEdge.right_face;
- END IF;
- ELSE -- ring on the right is a shell
- IF ST_Contains(ringRight, apoint) THEN
- RETURN closestEdge.right_face;
- ELSE
- RETURN closestEdge.left_face;
- END IF;
- END IF;
-
- RETURN 0; -- anything not in a face is in the universal face
-
-END;
-$BODY$ LANGUAGE 'plpgsql';
+RETURNS INT AS
+ 'MODULE_PATHNAME', 'GetFaceContainingPoint'
+LANGUAGE 'c' STABLE;
commit 6f9afe7e25c561ebcc0cdc009d305133cf2bb279
Author: Sandro Santilli <strk at kbt.io>
Date: Wed Aug 11 11:26:13 2021 +0200
Add ptarray_closest_segment_2d lwgeom API function
Includes unit test
diff --git a/liblwgeom/cunit/cu_ptarray.c b/liblwgeom/cunit/cu_ptarray.c
index 01b1f46..c393577 100644
--- a/liblwgeom/cunit/cu_ptarray.c
+++ b/liblwgeom/cunit/cu_ptarray.c
@@ -685,6 +685,43 @@ static void test_ptarray_closest_vertex_2d()
lwline_free(line);
}
+static void test_ptarray_closest_segment_2d()
+{
+ LWLINE *line;
+ POINTARRAY *pa;
+ double dist;
+ POINT2D qp;
+ const char *wkt;
+ int rv;
+
+ wkt = "LINESTRING (0 0 0, 1 0 0, 2 0 0, 3 0 10)";
+ line = lwgeom_as_lwline(lwgeom_from_text(wkt));
+ pa = line->points;
+
+ qp.x = qp.y = 0;
+ rv = ptarray_closest_segment_2d(pa, &qp, &dist);
+ ASSERT_INT_EQUAL(rv, 0);
+ ASSERT_DOUBLE_EQUAL(dist, 0);
+
+ qp.x = 1;
+ rv = ptarray_closest_segment_2d(pa, &qp, &dist);
+ ASSERT_INT_EQUAL(rv, 0);
+ ASSERT_DOUBLE_EQUAL(dist, 0);
+
+ qp.y = 1;
+ rv = ptarray_closest_segment_2d(pa, &qp, &dist);
+ ASSERT_INT_EQUAL(rv, 0);
+ ASSERT_DOUBLE_EQUAL(dist, 1);
+
+ qp.x = 5; qp.y = 0;
+ rv = ptarray_closest_segment_2d(pa, &qp, &dist);
+ ASSERT_INT_EQUAL(rv, 2);
+ ASSERT_DOUBLE_EQUAL(dist, 2);
+
+
+ lwline_free(line);
+}
+
/*
** Used by the test harness to register the tests in this file.
@@ -704,4 +741,5 @@ void ptarray_suite_setup(void)
PG_ADD_TEST(suite, test_ptarray_scale);
PG_ADD_TEST(suite, test_ptarray_scroll);
PG_ADD_TEST(suite, test_ptarray_closest_vertex_2d);
+ PG_ADD_TEST(suite, test_ptarray_closest_segment_2d);
}
diff --git a/liblwgeom/liblwgeom.h.in b/liblwgeom/liblwgeom.h.in
index f2a148c..c074d8b 100644
--- a/liblwgeom/liblwgeom.h.in
+++ b/liblwgeom/liblwgeom.h.in
@@ -1097,6 +1097,16 @@ extern POINTARRAY *ptarray_substring(POINTARRAY *pa, double d1, double d2,
extern int ptarray_closest_vertex_2d(const POINTARRAY *pa,
const POINT2D *qp, double *dist);
+/**
+ * @param pa the subject pointarray
+ * @param qp the query point
+ * @param dist optional output for actual distance from segment
+ * @return 0-based segment index for the closest segment
+ * (earliest segment in case of same distance)
+ */
+extern int ptarray_closest_segment_2d(const POINTARRAY *pa,
+ const POINT2D *qp, double *dist);
+
/**
diff --git a/liblwgeom/ptarray.c b/liblwgeom/ptarray.c
index db99148..b58c8b2 100644
--- a/liblwgeom/ptarray.c
+++ b/liblwgeom/ptarray.c
@@ -1296,6 +1296,39 @@ closest_point_on_segment(const POINT4D *p, const POINT4D *A, const POINT4D *B, P
}
int
+ptarray_closest_segment_2d(const POINTARRAY *pa, const POINT2D *qp, double *dist)
+{
+ const POINT2D *start = getPoint2d_cp(pa, 0), *end = NULL;
+ uint32_t t, seg=0;
+ double mindist=DBL_MAX;
+
+ /* Loop through pointarray looking for nearest segment */
+ for (t=1; t<pa->npoints; t++)
+ {
+ double dist_sqr;
+ end = getPoint2d_cp(pa, t);
+ dist_sqr = distance2d_sqr_pt_seg(qp, start, end);
+
+ if (dist_sqr < mindist)
+ {
+ mindist = dist_sqr;
+ seg=t-1;
+ if ( mindist == 0 )
+ {
+ LWDEBUG(3, "Breaking on mindist=0");
+ break;
+ }
+ }
+
+ start = end;
+ }
+
+ if ( dist ) *dist = sqrt(mindist);
+ return seg;
+}
+
+
+int
ptarray_closest_vertex_2d(const POINTARRAY *pa, const POINT2D *qp, double *dist)
{
uint32_t t, pn=0;
commit ec3f582abbf9d22d07338521317f8ab94d144c29
Author: Sandro Santilli <strk at kbt.io>
Date: Wed Aug 11 15:15:20 2021 +0200
Ensure getEdgeByNode correctly initializes the "geom" field
diff --git a/topology/postgis_topology.c b/topology/postgis_topology.c
index 7374e89..d23831d 100644
--- a/topology/postgis_topology.c
+++ b/topology/postgis_topology.c
@@ -809,6 +809,10 @@ fillEdgeFields(LWT_ISO_EDGE* edge, HeapTuple row, TupleDesc rowdesc, int fields)
edge->geom = NULL;
}
}
+ else
+ {
+ edge->geom = NULL;
+ }
}
static void
-----------------------------------------------------------------------
Summary of changes:
liblwgeom/cunit/cu_ptarray.c | 38 +++
liblwgeom/liblwgeom.h.in | 10 +
liblwgeom/liblwgeom_topo.h | 52 ++--
liblwgeom/liblwgeom_topo_internal.h | 4 +-
liblwgeom/lwgeom_topo.c | 261 +++++++++++++++++++--
liblwgeom/ptarray.c | 33 +++
topology/postgis_topology.c | 150 +++++++++---
topology/sql/query/GetFaceContainingPoint.sql.in | 287 +----------------------
8 files changed, 484 insertions(+), 351 deletions(-)
hooks/post-receive
--
PostGIS
More information about the postgis-tickets
mailing list