[SCM] PostGIS branch master updated. 3.6.0rc2-436-g0f7d3d562
git at osgeo.org
git at osgeo.org
Wed Apr 8 14:39:07 PDT 2026
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, master has been updated
via 0f7d3d5628d0cf598c27ffe69e6fb25e8d05fe5d (commit)
from 2c9277d3da0794f44dee60107d107cfd9939edeb (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 0f7d3d5628d0cf598c27ffe69e6fb25e8d05fe5d
Author: Paul Ramsey <pramsey at cleverelephant.ca>
Date: Wed Apr 8 14:36:49 2026 -0700
Clean up and add another special case to lw_dist2d_pt_seg()
When the seg is vertical or horizontal, we construct the
nearest points to avoid getting a difference between
the more-exact distance calculation that doesn't construct
and the naturally constructive process of building
a bbox, which can result in a situation where the
"real distance" between two objects can be less than
the "box distance", which is simply not allowed.
References #6026
diff --git a/liblwgeom/measures.c b/liblwgeom/measures.c
index cc3ea6800..e10d60bf9 100644
--- a/liblwgeom/measures.c
+++ b/liblwgeom/measures.c
@@ -2206,103 +2206,157 @@ End of New faster distance calculations
Functions in common for Brute force and new calculation
--------------------------------------------------------------------------------------------------------------*/
-/**
-lw_dist2d_comp from p to line A->B
-This one is now sending every occasion to lw_dist2d_pt_pt
-Before it was handling occasions where r was between 0 and 1 internally
-and just returning the distance without identifying the points.
-To get this points it was necessary to change and it also showed to be about 10%faster.
-*/
int
-lw_dist2d_pt_seg(const POINT2D *p, const POINT2D *A, const POINT2D *B, DISTPTS *dl)
+lw_dist2d_pt_seg(const POINT2D *C, const POINT2D *A, const POINT2D *B, DISTPTS *dl)
{
- double r;
-
- LWDEBUG(2, "lw_dist2d_pt_seg called");
-
- /*if start==end, then use pt distance */
- if ((A->x == B->x) && (A->y == B->y))
- {
- LWDEBUG(2, "lw_dist2d_pt_seg found first and last segment points being the same");
- return lw_dist2d_pt_pt(p, A, dl);
- }
-
+ POINT2D P;
+ int is_vertical = (A->x == B->x);
+ int is_horizontal = (A->y == B->y);
/*
- * otherwise, we use comp.graphics.algorithms
- * Frequently Asked Questions method
+ * If A == B, then use pt-pt distance
+ */
+ if (is_vertical && is_horizontal)
+ return lw_dist2d_pt_pt(C, A, dl);
+
+ /*
+ * We use comp.graphics.algorithms Frequently Asked Questions method
+ *
+ * Let the point be C (Cx,Cy) and the line be AB (Ax,Ay) to (Bx,By).
+ * Let P be the point of perpendicular projection of C on AB. The parameter
+ * r, which indicates P's position along AB, is computed by the dot product
+ * of AC and AB divided by the square of the length of AB:
+ *
+ * (1) AC dot AB
+ * r = ---------
+ * ||AB||^2
+ *
+ * The length of a line segment
+ *
+ * L = sqrt( (Bx-Ax)^2 + (By-Ay)^2 )
+ *
+ * So (1) expands to:
+ *
+ * (Cx-Ax)(Bx-Ax) + (Cy-Ay)(By-Ay)
+ * r = -------------------------------
+ * L^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
+ * 0<r<1 P is interior to AB
*
- * (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
- * 0<r<1 P is interior to AB
*/
- r = ((p->x - 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));
+ double dx = B->x - A->x;
+ double dy = B->y - A->y;
+ double L2 = (dx*dx) + (dy*dy);
- LWDEBUGF(2, "lw_dist2d_pt_seg found r = %.15g", r);
+ double r = ((C->x - A->x) * dx + (C->y - A->y) * dy) / L2;
- /*This is for finding the maxdistance.
- the maxdistance have to be between two vertices, compared to mindistance which can be between two vertices.*/
+ /*
+ * This is for finding the maxdistance.
+ * The maxdistance have to be between two vertices,
+ * compared to mindistance which can be between two vertices.
+ */
if (dl->mode == DIST_MAX)
{
if (r >= 0.5)
- return lw_dist2d_pt_pt(p, A, dl);
+ return lw_dist2d_pt_pt(C, A, dl);
else /* (r < 0.5) */
- return lw_dist2d_pt_pt(p, B, dl);
+ return lw_dist2d_pt_pt(C, B, dl);
}
- if (r < 0) /*If p projected on the line is outside point A*/
- return lw_dist2d_pt_pt(p, A, dl);
- if (r >= 1) /*If p projected on the line is outside point B or on point B*/
- return lw_dist2d_pt_pt(p, B, dl);
-
- /*If the point p is on the segment this is a more robust way to find out that*/
- if ((((A->y - p->y) * (B->x - A->x) == (A->x - p->x) * (B->y - A->y))) && (dl->mode == DIST_MIN))
- {
- lw_dist2d_distpts_set(dl, 0, p, p);
- }
-
-
/*
- (2)
- (Ay-Cy)(Bx-Ax)-(Ax-Cx)(By-Ay)
- s = -----------------------------
- L^2
+ * If projected point P is outside point A
+ */
+ if (r < 0)
+ return lw_dist2d_pt_pt(C, A, dl);
- Then the distance from C to P = |s|*L.
- */
+ /*
+ * If projected point P is is outside
+ * point B or on point B
+ */
+ if (r >= 1)
+ return lw_dist2d_pt_pt(C, B, dl);
- double s = ((A->y - p->y) * (B->x - A->x) - (A->x - p->x) * (B->y - A->y)) /
- ((B->x - A->x) * (B->x - A->x) + (B->y - A->y) * (B->y - A->y));
+ /*
+ * If point C is on the segment AB then
+ * distance is zero and nearest point is C
+ */
+ if ((A->y - C->y) * dx == (A->x - C->x) * dy)
+ {
+ lw_dist2d_distpts_set(dl, 0.0, C, C);
+ }
- double dist = fabs(s) * sqrt(((B->x - A->x) * (B->x - A->x) + (B->y - A->y) * (B->y - A->y)));
- if ( dist < dl->distance )
+ /*
+ * The point P can then be found:
+ *
+ * Px = Ax + r(Bx-Ax)
+ * Py = Ay + r(By-Ay)
+ *
+ * And the distance from A to P = r*L.
+ */
+
+ /*
+ * If the projection of point C on the segment is
+ * between A and B then we find that "point on segment"
+ * and send it to lw_dist2d_pt_pt
+ */
+ if (is_horizontal || is_vertical)
+ {
+ P.x = A->x + r * dx;
+ P.y = A->y + r * dy;
+ return lw_dist2d_pt_pt(C, &P, dl);
+ }
+
+ /*
+ *
+ * Use another parameter s to indicate the location along PC, with the
+ * following meaning:
+ * s<0 C is left of AB
+ * s>0 C is right of AB
+ * s=0 C is on AB
+ *
+ * Compute s as follows:
+ *
+ * (Ay-Cy)(Bx-Ax)-(Ax-Cx)(By-Ay)
+ * s = -----------------------------
+ * L^2
+ *
+ *
+ * Then the distance from C to P = s*L.
+ */
+
+ /*
+ * Calculate distance without reference to the
+ * projected point P.
+ */
+ double s = ((A->y - C->y) * dx - (A->x - C->x) * dy) / L2;
+
+ double dist = fabs(s) * sqrt(L2);
+ // double dist = fabs(s) * hypot(dx, dy);
+ if (dist < dl->distance)
{
dl->distance = dist;
{
- POINT2D c;
- c.x = A->x + r * (B->x - A->x);
- c.y = A->y + r * (B->y - A->y);
+ P.x = A->x + r * dx;
+ P.y = A->y + r * dy;
if (dl->twisted > 0)
{
- dl->p1 = *p;
- dl->p2 = c;
+ dl->p1 = *C;
+ dl->p2 = P;
}
else
{
- dl->p1 = c;
- dl->p2 = *p;
+ dl->p1 = P;
+ dl->p2 = *C;
}
}
}
-
return LW_TRUE;
}
diff --git a/regress/core/knn_recheck.sql b/regress/core/knn_recheck.sql
index f08b7f0d6..45278b7fc 100644
--- a/regress/core/knn_recheck.sql
+++ b/regress/core/knn_recheck.sql
@@ -239,3 +239,31 @@ SET enable_seqscan = false;
SELECT '#5782', l FROM t5782 ORDER BY g <-> 'POINT(18.006691126034692 69.04048768506776)'::geometry;
DROP TABLE t5782;
SET enable_seqscan to default;
+
+--
+-- https://trac.osgeo.org/postgis/ticket/6026
+--
+CREATE TABLE reproduce_6026 (
+ id SERIAL PRIMARY KEY,
+ geometry GEOMETRY(POLYGON, 28992)
+);
+INSERT INTO reproduce_6026 (geometry) VALUES
+ (ST_GeomFromText('POLYGON ((133894.531 490000, 132700 490000, 132460 490000, 131703.958 489011.71, 131695 489000, 132880 489000, 132880 488416.725, 133206.605 488583.011, 134029.81 489022.36, 133917.695 490000, 133894.531 490000))', 28992)),
+ (ST_GeomFromText('POLYGON ((138000 492000, 141276.49 492000, 141050 492550, 140967.11 492751.305, 139951.651 495217.42, 139490.072 495217.42, 136358.892 495217.42, 135900 495217.42, 135301.1 495217.42, 135118.887 495296.573, 135115.873 495264.338, 135106.301 495164.456, 135098.485 495065.802, 135090.476 494967.014, 135085.764 494917.296, 135078.439 494866.459, 135070.381 494817.483, 135060.353 494767.172, 135036.351 494674.319, 135017.779 494615.505, 134991.438 494544.049, 134972.66 494501.055, 134948.148 494446.666, 134922.544 494398.914, 134872.859 494311.273, 134859.929 494292.444, 134843.402 494264.587, 134815.107 494200.454, 134770.855 494120.756, 134728.422 494053.048, 134932.74 493923.77, 133708.3 491825.92, 133917.695 490000, 134029.81 489022.36, 138000 492000))', 28992)),
+ (ST_GeomFromText('POLYGON ((131215.201 487105.005, 134029.81 486992.15, 134029.81 489022.36, 133206.605 488583.011, 132880 488416.725, 132198.685 488069.845, 131375.721 487651.916, 131199.012 487560.733, 131198.76 487560.603, 131215.201 487105.005))', 28992)),
+ (ST_GeomFromText('POLYGON ((134029.81 489022.36, 136889.1 488581.7, 138000 492000, 134029.81 489022.36))', 28992)),
+ (ST_GeomFromText('POLYGON ((133917.695 490000, 133708.3 491825.92, 133541.312 491951.805, 133460 491840, 132820 491000, 132210 490190, 132460 490000, 132700 490000, 133894.531 490000, 133917.695 490000))', 28992));
+
+-- create index
+CREATE INDEX reproduce_6026_x
+ON reproduce_6026 USING GIST(geometry);
+-- force index usage
+SET enable_seqscan = OFF;
+-- Reproduce error
+SELECT '#6026', id
+ FROM reproduce_6026
+ ORDER BY geometry <-> ST_SetSRID(ST_MakePoint(133718, 489222), 28992)
+ LIMIT 4;
+
+DROP TABLE reproduce_6026;
+
diff --git a/regress/core/knn_recheck_expected b/regress/core/knn_recheck_expected
index 226f99ed7..051131700 100644
--- a/regress/core/knn_recheck_expected
+++ b/regress/core/knn_recheck_expected
@@ -107,3 +107,7 @@
#3418|0.5500000|0.5500000
#5782|A
#5782|B
+#6026|1
+#6026|2
+#6026|3
+#6026|4
-----------------------------------------------------------------------
Summary of changes:
liblwgeom/measures.c | 176 +++++++++++++++++++++++++-------------
regress/core/knn_recheck.sql | 28 ++++++
regress/core/knn_recheck_expected | 4 +
3 files changed, 147 insertions(+), 61 deletions(-)
hooks/post-receive
--
PostGIS
More information about the postgis-tickets
mailing list