[postgis-tickets] [SCM] PostGIS branch main updated. 3.1.0rc1-397-g6571dc8
git at osgeo.org
git at osgeo.org
Thu Aug 12 14:06:53 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 6571dc8fde5cb396bfaf6db394f39b765535028c (commit)
from b4f0beefdc3bcd2ed555ce9da6f096655e875ec3 (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 6571dc8fde5cb396bfaf6db394f39b765535028c
Author: Sandro Santilli <strk at kbt.io>
Date: Thu Aug 12 21:13:25 2021 +0200
Use angles to determine side of point on edge when closest is a vertex
Closes #4962 again
Includes regress test for getFaceContainingPoint
(failing without the change)
diff --git a/liblwgeom/lwgeom_topo.c b/liblwgeom/lwgeom_topo.c
index 23b99c1..49aa52b 100644
--- a/liblwgeom/lwgeom_topo.c
+++ b/liblwgeom/lwgeom_topo.c
@@ -7182,6 +7182,8 @@ lwt_GetFaceContainingPoint(LWT_TOPOLOGY* topo, const LWPOINT* pt)
lwerror("error computing azimuth of shortestLine [%.15g %.15g,%.15g %.15g]",
shortestLineP0->x, shortestLineP0->y,
shortestLineP1->x, shortestLineP1->y);
+ lwgeom_free(shortestLine);
+ _lwt_release_edges(closestEdge, 1);
return -1;
}
@@ -7193,10 +7195,12 @@ lwt_GetFaceContainingPoint(LWT_TOPOLOGY* topo, const LWPOINT* pt)
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);
+ lwgeom_free(shortestLine);
_lwt_release_edges(closestEdge, 1);
return -1;
}
+ lwgeom_free(shortestLine);
_lwt_release_edges(closestEdge, 1);
return ee.cwFace;
@@ -7209,6 +7213,7 @@ lwt_GetFaceContainingPoint(LWT_TOPOLOGY* topo, const LWPOINT* pt)
if ( closestEdge->face_left == closestEdge->face_right )
{
containingFace = closestEdge->face_left;
+ lwgeom_free(shortestLine);
_lwt_release_edges(closestEdge, 1);
return containingFace;
}
@@ -7217,6 +7222,7 @@ lwt_GetFaceContainingPoint(LWT_TOPOLOGY* topo, const LWPOINT* pt)
if ( dist == 0 )
{
/* We checked the dangling case above */
+ lwgeom_free(shortestLine);
_lwt_release_edges(closestEdge, 1);
lwerror("Two or more faces found");
return -1;
@@ -7225,11 +7231,11 @@ lwt_GetFaceContainingPoint(LWT_TOPOLOGY* topo, const LWPOINT* pt)
/* 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);
+ shortestLineP1 = getPoint2d_cp(((LWLINE *)shortestLine)->points, 1);
+ closestSegmentIndex = ptarray_closest_segment_2d(closestEdge->geom->points, shortestLineP1, 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);
+ closestSegmentP1 = getPoint2d_cp(closestEdge->geom->points, closestSegmentIndex + 1);
LWDEBUGF(1, "Closest segment to edge %d is LINESTRING(%g %g, %g %g)",
closestEdge->edge_id,
closestSegmentP0->x,
@@ -7239,11 +7245,94 @@ lwt_GetFaceContainingPoint(LWT_TOPOLOGY* topo, const LWPOINT* pt)
);
/* Find on which side of the segment the query point lays */
+ shortestLineP0 = getPoint2d_cp(((LWLINE *)shortestLine)->points, 0);
+ if ( p2d_same(shortestLineP0, closestSegmentP0) )
+ {
+ /* Closest point is first point of closest segment (this should
+ * never happen as we'd have returned the previous segment in
+ * this case, unless it was the very first segment, and in that
+ * case we'd have cought this earlier, as it would be a NODE
+ *
+ * This would be more of an assert...
+ */
+ lwerror("Closest point is first point of closest segment, unexpectedly");
+ lwgeom_free(shortestLine);
+ _lwt_release_edges(closestEdge, 1);
+ return -1;
+ }
+ else if ( p2d_same(shortestLineP0, closestSegmentP1) )
+ {
+ /* Closest point is last point of the closest segment,
+ * so we need to check if rotating the closest segment
+ * clockwise around its last point we encounter
+ * the shortestLine first (which means it's on the left
+ * of the closest edge) or the next segment first (which
+ * means the shortestLine is on the right)
+ */
+ const POINT2D *nextSegmentP0, *nextSegmentP1;
+ double azS0, azS1, azSL;
- do {
- containingFace = -1;
+ nextSegmentP0 = closestSegmentP1;
+ /* This would be more of an assert...*/
+ if ( closestSegmentIndex + 1 > closestEdge->geom->points->npoints )
+ {
+ lwerror("closestSegmentIndex is unexpectedly the last one and we didn't exit earlier as it would be a node");
+ return -1;
+ }
+ nextSegmentP1 = getPoint2d_cp(closestEdge->geom->points, closestSegmentIndex + 2);
+
+ if ( ! azimuth_pt_pt(closestSegmentP1, closestSegmentP0, &azS0)) {
+ lwerror("error computing azimuth of reversse closest segment [%.15g %.15g,%.15g %.15g]",
+ closestSegmentP1->x, closestSegmentP1->y,
+ closestSegmentP0->x, closestSegmentP0->y);
+ lwgeom_free(shortestLine);
+ _lwt_release_edges(closestEdge, 1);
+ return -1;
+ }
+ if ( ! azimuth_pt_pt(nextSegmentP0, nextSegmentP1, &azS1)) {
+ lwerror("error computing azimuth of next segment [%.15g %.15g,%.15g %.15g]",
+ nextSegmentP0->x, nextSegmentP0->y,
+ nextSegmentP1->x, nextSegmentP1->y);
+ lwgeom_free(shortestLine);
+ _lwt_release_edges(closestEdge, 1);
+ return -1;
+ }
+ if ( ! azimuth_pt_pt(shortestLineP0, shortestLineP1, &azSL) ) {
+ lwerror("error computing azimuth of shortestLine [%.15g %.15g,%.15g %.15g]",
+ shortestLineP0->x, shortestLineP0->y,
+ shortestLineP1->x, shortestLineP1->y);
+ lwgeom_free(shortestLine);
+ _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 last point of closest segment: S0:%g, S1:%g (+%g), SL:%g (+%g)",
+ azS0,
+ azS1, angle_S0_S1,
+ azSL, angle_S0_SL
+ );
+ if ( angle_S0_SL < angle_S0_S1 )
+ {
+ /* shortest line was encountered first, is on the left */
+ containingFace = closestEdge->face_left;
+ }
+ else
+ {
+ /* shortest line 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 */
- shortestLineP1 = getPoint2d_cp(((LWLINE *)shortestLine)->points, 1);
LWDEBUGF(1, "Calling lw_segment_side((%g,%g),(%g,%g),(%g,%g)",
closestSegmentP0->x,
closestSegmentP0->y,
@@ -7259,26 +7348,20 @@ lwt_GetFaceContainingPoint(LWT_TOPOLOGY* topo, const LWPOINT* pt)
if ( closestSegmentSide == -1 ) /* left */
{
containingFace = closestEdge->face_left;
- break;
}
else if ( closestSegmentSide == 1 ) /* right */
{
containingFace = closestEdge->face_right;
- break;
+ }
+ else
+ {
+ lwerror("Unexpected collinearity reported from lw_segment_side");
+ _lwt_release_edges(closestEdge, 1);
+ lwgeom_free(shortestLine);
+ return -1;
}
- /* 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);
diff --git a/topology/test/regress/getfacecontainingpoint.sql b/topology/test/regress/getfacecontainingpoint.sql
index 9763a2d..d4befbd 100644
--- a/topology/test/regress/getfacecontainingpoint.sql
+++ b/topology/test/regress/getfacecontainingpoint.sql
@@ -51,6 +51,18 @@ SELECT 't8', topology.GetFaceContainingPoint('city_data', 'POINT(24 16.9)');
-- Query point collinear with closest incoming segment (right of next one)
SELECT 't9', topology.GetFaceContainingPoint('city_data', 'POINT(7.1 31)');
+-- Query point on the right of the closest segment but left of ring
+SELECT 't10', topology.GetFaceContainingPoint('city_data', 'POINT(26 16.8)');
+
+-- Query point on the left of the closest segment and left of ring
+SELECT 't11', topology.GetFaceContainingPoint('city_data', 'POINT(26 17.2)');
+
+-- Query point on the left of the closest segment but right of ring
+SELECT 't12', topology.GetFaceContainingPoint('city_data', 'POINT(33.1 17.8)');
+
+-- Query point on the right of the closest segment and right of ring
+SELECT 't13', topology.GetFaceContainingPoint('city_data', 'POINT(33.1 17.2)');
+
-- Query point on a node incident to multiple faces
SELECT 'e1', topology.GetFaceContainingPoint('city_data', 'POINT(21 14)');
diff --git a/topology/test/regress/getfacecontainingpoint_expected b/topology/test/regress/getfacecontainingpoint_expected
index 251363c..ef034b2 100644
--- a/topology/test/regress/getfacecontainingpoint_expected
+++ b/topology/test/regress/getfacecontainingpoint_expected
@@ -18,6 +18,10 @@ t6|0
t7|11
t8|4
t9|1
+t10|10
+t11|10
+t12|10
+t13|10
ERROR: Two or more faces found
ERROR: Two or more faces found
ERROR: Two or more faces found
-----------------------------------------------------------------------
Summary of changes:
liblwgeom/lwgeom_topo.c | 123 +++++++++++++++++----
topology/test/regress/getfacecontainingpoint.sql | 12 ++
.../test/regress/getfacecontainingpoint_expected | 4 +
3 files changed, 119 insertions(+), 20 deletions(-)
hooks/post-receive
--
PostGIS
More information about the postgis-tickets
mailing list