[SCM] PostGIS branch master updated. 3.7.0rc1-103-g83e528504
git at osgeo.org
git at osgeo.org
Fri Oct 9 01:35:02 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 83e52850466d7566ae88f0ee832829cbd71a3ce1 (commit)
from 585ba4895376d5300bf9ec08ea5aef1fd8a2b728 (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 83e52850466d7566ae88f0ee832829cbd71a3ce1
Author: Sandro Santilli <strk at kbt.io>
Date: Fri Oct 9 07:44:29 2026 +0000
Use exact arithmetic for edge crossing checks, move it to topology
Fixes #6143, including unit and regression test for it
diff --git a/liblwgeom/cunit/cu_algorithm.c b/liblwgeom/cunit/cu_algorithm.c
index 9c2434e4a..e98b571c2 100644
--- a/liblwgeom/cunit/cu_algorithm.c
+++ b/liblwgeom/cunit/cu_algorithm.c
@@ -1880,158 +1880,6 @@ static void test_ptarray_simplify(void)
lwgeom_free(geom2);
}
-/*
-** Assert the intersection status of a pair of lines described by WKT.
-**
-** The check must be symmetric, so both argument orders are verified.
-*/
-#define ASSERT_NONBOUNDARY_INTERSECTION(wkt1, wkt2, expected) \
- do \
- { \
- LWLINE *l1 = (LWLINE *)lwgeom_from_wkt(wkt1, LW_PARSER_CHECK_NONE); \
- LWLINE *l2 = (LWLINE *)lwgeom_from_wkt(wkt2, LW_PARSER_CHECK_NONE); \
- int exp = (expected); \
- CU_ASSERT_EQUAL(lwline_have_nonboundary_2d_intersection(l1, l2), exp); \
- CU_ASSERT_EQUAL(lwline_have_nonboundary_2d_intersection(l2, l1), exp); \
- lwline_free(l1); \
- lwline_free(l2); \
- } while (0)
-
-static void test_lwline_have_nonboundary_2d_intersection(void)
-{
- /* Non intersecting lines */
- ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 100 0)", "LINESTRING(200 0, 400 0)", 0);
- ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 100 0)", "LINESTRING(0 50, 100 50)", 0);
- /* The two lines would cross, but only outside of both segments */
- ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 100 0)", "LINESTRING(200 100, 200 -100)", 0);
- /* The two lines would touch, but only outside of both segments */
- ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 100 0)", "LINESTRING(-100 0, -50 0)", 0);
- /* Lines parallel and apart, also on a diagonal */
- ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 100 100)", "LINESTRING(0 50, 100 150)", 0);
-
- /* Touching on a boundary point of both lines */
- /* End to start */
- ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 100 0)", "LINESTRING(100 0, 200 0)", 0);
- /* Start to start */
- ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 100 0)", "LINESTRING(0 0, 0 100)", 0);
- /* Start to end */
- ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 100 0)", "LINESTRING(100 0, 100 100)", 0);
- /* End to end */
- ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 100 0)", "LINESTRING(200 0, 100 0)", 0);
- /* Same point, opposite directions */
- ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(100 0, 0 0)", "LINESTRING(0 0, -100 0)", 0);
- /* Collinear, touching end to start on a vertical line */
- ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 0 100)", "LINESTRING(0 100, 0 200)", 0);
- /* Collinear, touching end to start on a diagonal */
- ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 100 100)", "LINESTRING(100 100, 200 200)", 0);
- /* Collinear, touching on the interior vertex of the first line */
- ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 50 0, 100 0)", "LINESTRING(100 0, 150 0)", 0);
- /* Diverging from the shared start point */
- ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 100 0)", "LINESTRING(0 0, 50 50)", 0);
- /* The endpoint of the second line is the start point of the first */
- ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 100 0)", "LINESTRING(-50 0, 0 0)", 0);
- /* Zero length second line, on the boundary of the first one */
- ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 100 0)", "LINESTRING(100 0, 100 0)", 0);
- /* Zero length lines, on the same point */
- ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(50 50, 50 50)", "LINESTRING(50 50, 50 50)", 0);
- /* Closed line, touching another line on its boundary point */
- ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 100 0, 100 100, 0 100, 0 0)",
- "LINESTRING(0 0, -100 0)",
- 0);
-
- /* Crossing lines */
- ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 100 0)", "LINESTRING(50 10, 50 -10)", 1);
- /* Crossing on a diagonal */
- ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 100 100)", "LINESTRING(0 100, 100 0)", 1);
- /* Crossing, but one line is zero length and on the other one interior */
- ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 100 0)", "LINESTRING(50 0, 50 0)", 1);
-
- /* Identical lines */
- ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 100 0)", "LINESTRING(0 0, 100 0)", 1);
- ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 50 0, 100 0)", "LINESTRING(0 0, 50 0, 100 0)", 1);
- /* Identical lines, opposite direction, on a vertical */
- ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 0 100)", "LINESTRING(0 100, 0 0)", 1);
-
- /* Endpoint of one line on the interior of the other one */
- ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 100 0)", "LINESTRING(50 0, 50 50)", 1);
- ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 100 0)", "LINESTRING(50 -50, 50 0)", 1);
- /* Zero length first line, on the interior of the second one */
- ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(50 0, 50 0)", "LINESTRING(0 0, 100 0)", 1);
- /* Collinear, zero length second line on the interior of the first one */
- ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 50 0, 100 0)", "LINESTRING(50 0, 50 0)", 1);
- /* Interior vertex of the first line, endpoint of the second one */
- ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 50 0, 100 0)", "LINESTRING(50 0, 50 100)", 1);
- /* Interior vertex of the first line, interior of the second one */
- ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 50 0, 100 0)", "LINESTRING(50 -50, 50 50)", 1);
- /* Collinear, overlapping for a positive length */
- ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 100 0)", "LINESTRING(50 0, 150 0)", 1);
- /* Collinear, overlapping for a positive length, opposite direction */
- ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(100 0, 0 0)", "LINESTRING(-50 0, 50 0)", 1);
- /* Collinear, second line contained in the first one */
- ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 100 0)", "LINESTRING(25 0, 75 0)", 1);
- /* Collinear, overlapping on a diagonal */
- ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 100 100)", "LINESTRING(50 50, 150 150)", 1);
- /* Collinear, overlapping on a vertical */
- ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 0 100)", "LINESTRING(0 50, 0 200)", 1);
- /* Collinear, overlapping, the shared point is on the line boundary */
- ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 100 0)", "LINESTRING(0 0, 50 0)", 1);
- /* Collinear, overlapping across multiple segments */
- ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 100 0)", "LINESTRING(-50 0, 25 0, 150 0)", 1);
- /* Collinear, the second line is back and forth over the first one */
- ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 100 0)", "LINESTRING(50 0, 100 100, 50 0)", 1);
- /* Second line touches an interior vertex of first line */
- ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 100 0, 200 0)", "LINESTRING(100 0, 100 -50)", 1);
-
- /*
- ** Multiple intersections on the same pair of lines, the shared
- ** boundary point must not hide the non-boundary ones
- */
- ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 100 0)", "LINESTRING(0 0, 50 50, 50 -50, 0 0)", 1);
- ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 100 0)",
- "LINESTRING(0 0, 50 50, 100 0, 50 -50, 0 0)",
- 1);
-
- /* Z and M coordinates are ignored */
- ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING Z (0 0 100, 100 0 200)",
- "LINESTRING M (50 10 1, 50 -10 2)",
- 1);
- ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING Z (0 0 1, 100 0 1)", "LINESTRING Z (0 0 2, 100 0 2)", 1);
- ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING ZM (0 0 1 1, 100 0 1 1)",
- "LINESTRING Z (100 0 9, 100 0 9)",
- 0);
-
- /* Empty lines have nothing to intersect */
- ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING EMPTY", "LINESTRING(0 0, 100 0)", 0);
- ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 100 0)", "LINESTRING EMPTY", 0);
- ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING EMPTY", "LINESTRING EMPTY", 0);
-
- /* Single point lines have no segment to intersect */
- {
- POINT4D pt;
- LWLINE *single, *line;
-
- pt.x = 0;
- pt.y = 0;
- pt.z = 0;
- pt.m = 0;
-
- single = lwline_construct_empty(SRID_UNKNOWN, 0, 0);
- ptarray_append_point(single->points, &pt, LW_TRUE);
- line = (LWLINE*)lwgeom_from_wkt("LINESTRING(0 0, 100 0)", LW_PARSER_CHECK_NONE);
-
- CU_ASSERT_EQUAL( lwline_have_nonboundary_2d_intersection(single, line), 0 );
- CU_ASSERT_EQUAL( lwline_have_nonboundary_2d_intersection(line, single), 0 );
- CU_ASSERT_EQUAL( lwline_have_nonboundary_2d_intersection(single, single), 0 );
-
- lwline_free(single);
- lwline_free(line);
- }
-
- /* NULL lines are handled gracefully */
- CU_ASSERT_EQUAL( lwline_have_nonboundary_2d_intersection(NULL, NULL), 0 );
-}
-
-
/*
** Used by test harness to register the tests in this file.
*/
@@ -2071,5 +1919,4 @@ void algorithms_suite_setup(void)
PG_ADD_TEST(suite,test_lwpoly_construct_circle);
PG_ADD_TEST(suite,test_trim_bits);
PG_ADD_TEST(suite,test_lwgeom_remove_repeated_points);
- PG_ADD_TEST(suite,test_lwline_have_nonboundary_2d_intersection);
}
diff --git a/liblwgeom/cunit/cu_topo.c b/liblwgeom/cunit/cu_topo.c
index 490b95058..41dfe9cdf 100644
--- a/liblwgeom/cunit/cu_topo.c
+++ b/liblwgeom/cunit/cu_topo.c
@@ -13,6 +13,7 @@
#include "liblwgeom_internal.h"
#include "topo/liblwgeom_topo.h"
+#include "lwgeom_geos.h"
#include "cu_tester.h"
static LWGEOM *
@@ -47,9 +48,169 @@ test_lwt_IsTopoRingCCW_large_finite_coordinates(void)
lwgeom_free(geom);
}
+/*
+** Assert the non-boundary intersection status of a pair of lines
+** described by WKT, as reported by the exact arithmetic predicate
+** used by the topology edge crossing checks.
+**
+** The check must be symmetric, so both argument orders are verified.
+*/
+#define ASSERT_NONBOUNDARY_INTERSECTION(wkt1, wkt2, expected) \
+ do \
+ { \
+ LWLINE *l1 = (LWLINE *)lwgeom_from_text(wkt1); \
+ LWLINE *l2 = (LWLINE *)lwgeom_from_text(wkt2); \
+ int exp = (expected); \
+ CU_ASSERT_EQUAL(lwt_LineHaveNonBoundary2DIntersection(l1, l2), exp); \
+ CU_ASSERT_EQUAL(lwt_LineHaveNonBoundary2DIntersection(l2, l1), exp); \
+ lwline_free(l1); \
+ lwline_free(l2); \
+ } while (0)
+
+static void
+test_lwt_LineHaveNonBoundary2DIntersection(void)
+{
+ /* Non intersecting lines */
+ ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 100 0)", "LINESTRING(200 0, 400 0)", 0);
+ ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 100 0)", "LINESTRING(0 50, 100 50)", 0);
+ /* The two lines would cross, but only outside of both segments */
+ ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 100 0)", "LINESTRING(200 100, 200 -100)", 0);
+ /* The two lines would touch, but only outside of both segments */
+ ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 100 0)", "LINESTRING(-100 0, -50 0)", 0);
+ /* Lines parallel and apart, also on a diagonal */
+ ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 100 100)", "LINESTRING(0 50, 100 150)", 0);
+
+ /* Touching on a boundary point of both lines */
+ /* End to start */
+ ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 100 0)", "LINESTRING(100 0, 200 0)", 0);
+ /* Start to start */
+ ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 100 0)", "LINESTRING(0 0, 0 100)", 0);
+ /* Start to end */
+ ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 100 0)", "LINESTRING(100 0, 100 100)", 0);
+ /* End to end */
+ ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 100 0)", "LINESTRING(200 0, 100 0)", 0);
+ /* Same point, opposite directions */
+ ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(100 0, 0 0)", "LINESTRING(0 0, -100 0)", 0);
+ /* Collinear, touching end to start on a vertical line */
+ ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 0 100)", "LINESTRING(0 100, 0 200)", 0);
+ /* Collinear, touching end to start on a diagonal */
+ ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 100 100)", "LINESTRING(100 100, 200 200)", 0);
+ /* Collinear, touching on the interior vertex of the first line */
+ ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 50 0, 100 0)", "LINESTRING(100 0, 150 0)", 0);
+ /* Diverging from the shared start point */
+ ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 100 0)", "LINESTRING(0 0, 50 50)", 0);
+ /* The endpoint of the second line is the start point of the first */
+ ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 100 0)", "LINESTRING(-50 0, 0 0)", 0);
+ /* Zero length second line, on the boundary of the first one */
+ ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 100 0)", "LINESTRING(100 0, 100 0)", 0);
+ /* Zero length lines, on the same point */
+ ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(50 50, 50 50)", "LINESTRING(50 50, 50 50)", 0);
+ /* Closed line, touching another line on its boundary point */
+ ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 100 0, 100 100, 0 100, 0 0)",
+ "LINESTRING(0 0, -100 0)",
+ 0);
+
+ /* Crossing lines */
+ ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 100 0)", "LINESTRING(50 10, 50 -10)", 1);
+ /* Crossing on a diagonal */
+ ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 100 100)", "LINESTRING(0 100, 100 0)", 1);
+ /* Crossing, but one line is zero length and on the other one interior */
+ ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 100 0)", "LINESTRING(50 0, 50 0)", 1);
+
+ /* Identical lines */
+ ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 100 0)", "LINESTRING(0 0, 100 0)", 1);
+ ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 50 0, 100 0)", "LINESTRING(0 0, 50 0, 100 0)", 1);
+ /* Identical lines, opposite direction, on a vertical */
+ ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 0 100)", "LINESTRING(0 100, 0 0)", 1);
+
+ /* Endpoint of one line on the interior of the other one */
+ ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 100 0)", "LINESTRING(50 0, 50 50)", 1);
+ ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 100 0)", "LINESTRING(50 -50, 50 0)", 1);
+ /* Zero length first line, on the interior of the second one */
+ ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(50 0, 50 0)", "LINESTRING(0 0, 100 0)", 1);
+ /* Collinear, zero length second line on the interior of the first one */
+ ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 50 0, 100 0)", "LINESTRING(50 0, 50 0)", 1);
+ /* Interior vertex of the first line, endpoint of the second one */
+ ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 50 0, 100 0)", "LINESTRING(50 0, 50 100)", 1);
+ /* Interior vertex of the first line, interior of the second one */
+ ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 50 0, 100 0)", "LINESTRING(50 -50, 50 50)", 1);
+ /* Collinear, overlapping for a positive length */
+ ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 100 0)", "LINESTRING(50 0, 150 0)", 1);
+ /* Collinear, overlapping for a positive length, opposite direction */
+ ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(100 0, 0 0)", "LINESTRING(-50 0, 50 0)", 1);
+ /* Collinear, second line contained in the first one */
+ ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 100 0)", "LINESTRING(25 0, 75 0)", 1);
+ /* Collinear, overlapping on a diagonal */
+ ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 100 100)", "LINESTRING(50 50, 150 150)", 1);
+ /* Collinear, overlapping on a vertical */
+ ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 0 100)", "LINESTRING(0 50, 0 200)", 1);
+ /* Collinear, overlapping, the shared point is on the line boundary */
+ ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 100 0)", "LINESTRING(0 0, 50 0)", 1);
+ /* Collinear, overlapping across multiple segments */
+ ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 100 0)", "LINESTRING(-50 0, 25 0, 150 0)", 1);
+ /* Collinear, the second line is back and forth over the first one */
+ ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 100 0)", "LINESTRING(50 0, 100 100, 50 0)", 1);
+ /* Second line touches an interior vertex of first line */
+ ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 100 0, 200 0)", "LINESTRING(100 0, 100 -50)", 1);
+
+ /*
+ ** Multiple intersections on the same pair of lines, the shared
+ ** boundary point must not hide the non-boundary ones
+ */
+ ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 100 0)", "LINESTRING(0 0, 50 50, 50 -50, 0 0)", 1);
+ ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 100 0)",
+ "LINESTRING(0 0, 50 50, 100 0, 50 -50, 0 0)",
+ 1);
+
+ /* Z and M coordinates are ignored */
+ ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING Z (0 0 100, 100 0 200)",
+ "LINESTRING M (50 10 1, 50 -10 2)",
+ 1);
+ ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING Z (0 0 1, 100 0 1)", "LINESTRING Z (0 0 2, 100 0 2)", 1);
+ ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING ZM (0 0 1 1, 100 0 1 1)",
+ "LINESTRING Z (100 0 9, 100 0 9)",
+ 0);
+
+ /* Empty lines have nothing to intersect */
+ ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING EMPTY", "LINESTRING(0 0, 100 0)", 0);
+ ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(0 0, 100 0)", "LINESTRING EMPTY", 0);
+ ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING EMPTY", "LINESTRING EMPTY", 0);
+
+ /* Single point lines have no segment to intersect */
+ {
+ POINT4D pt;
+ LWLINE *single, *line;
+
+ pt.x = 0;
+ pt.y = 0;
+ pt.z = 0;
+ pt.m = 0;
+
+ single = lwline_construct_empty(SRID_UNKNOWN, 0, 0);
+ ptarray_append_point(single->points, &pt, LW_TRUE);
+ line = (LWLINE *)lwgeom_from_text("LINESTRING(0 0, 100 0)");
+
+ CU_ASSERT_EQUAL( lwt_LineHaveNonBoundary2DIntersection(single, line), 0 );
+ CU_ASSERT_EQUAL( lwt_LineHaveNonBoundary2DIntersection(line, single), 0 );
+ CU_ASSERT_EQUAL( lwt_LineHaveNonBoundary2DIntersection(single, single), 0 );
+
+ lwline_free(single);
+ lwline_free(line);
+ }
+
+ /* NULL lines are handled gracefully */
+ CU_ASSERT_EQUAL( lwt_LineHaveNonBoundary2DIntersection(NULL, NULL), 0 );
+
+ /* The pair of lines below is disjoint
+ * See https://trac.osgeo.org/postgis/ticket/6143 */
+ ASSERT_NONBOUNDARY_INTERSECTION("LINESTRING(-1e12 -1e12, -1 1e-5)", "LINESTRING(1e5 -1e5, -1 -1e-11)", 0);
+
+}
+
void
topo_suite_setup(void)
{
CU_pSuite suite = CU_add_suite("topology", NULL, NULL);
PG_ADD_TEST(suite, test_lwt_IsTopoRingCCW_large_finite_coordinates);
+ PG_ADD_TEST(suite, test_lwt_LineHaveNonBoundary2DIntersection);
}
diff --git a/liblwgeom/liblwgeom_internal.h b/liblwgeom/liblwgeom_internal.h
index 179ac4460..0e8cd9384 100644
--- a/liblwgeom/liblwgeom_internal.h
+++ b/liblwgeom/liblwgeom_internal.h
@@ -481,7 +481,23 @@ int ptarray_raycast_intersections(const POINTARRAY *pa, const POINT2D *p, int *o
int ptarrayarc_raycast_intersections(const POINTARRAY *pa, const POINT2D *p, int *on_boundary);
int lwcompound_contains_point(const LWCOMPOUND *comp, const POINT2D *pt);
int lwgeom_contains_point(const LWGEOM *geom, const POINT2D *pt);
-int lwline_have_nonboundary_2d_intersection(const LWLINE *line1, const LWLINE *line2);
+
+/*
+* The ways two line segments can intersect.
+* Returned by the exact segment intersection test used by the
+* topology module (lwt_segment_intersection_exact in lwgeom_topo.c).
+*/
+enum LWSEGMENT_INTERSECTION
+{
+ /* The segments do not intersect */
+ LWSEG_INTERSECT_NONE = 0,
+ /* The segments are collinear and overlap for a positive length */
+ LWSEG_INTERSECT_OVERLAP,
+ /* The segments cross at a point interior to both of them */
+ LWSEG_INTERSECT_CROSS,
+ /* The segments touch at a single point, an endpoint of one of them */
+ LWSEG_INTERSECT_POINT
+};
/**
* Split a line by a point and push components to the provided multiline.
diff --git a/liblwgeom/lwline.c b/liblwgeom/lwline.c
index d8398dfc0..50185f510 100644
--- a/liblwgeom/lwline.c
+++ b/liblwgeom/lwline.c
@@ -747,264 +747,3 @@ lwline_extend(const LWLINE *line, double distance_forward, double distance_backw
}
return lwline_construct(line->srid, NULL, opa);
}
-
-/*
-* The ways two line segments can intersect.
-* Returned by lwsegment_intersection()
-*/
-enum LWSEGMENT_INTERSECTION
-{
- /* The segments do not intersect */
- LWSEG_INTERSECT_NONE = 0,
- /* The segments are collinear and overlap for a positive length */
- LWSEG_INTERSECT_OVERLAP,
- /* The segments cross at a point interior to both of them */
- LWSEG_INTERSECT_CROSS,
- /* The segments touch at a single point, an endpoint of one of them */
- LWSEG_INTERSECT_POINT
-};
-
-/**
-* Does the point lie on the segment, endpoints included?
-*
-* Unlike lw_pt_on_segment(), which excludes the segment endpoints.
-*
-* @param A1 first point of the segment
-* @param A2 second point of the segment
-* @param P the point to test
-*
-* @return LW_TRUE if P is on the A1/A2 segment
-*/
-static int
-lw_pt_on_segment_inclusive(const POINT2D *A1, const POINT2D *A2, const POINT2D *P)
-{
- if ( lw_segment_side(A1, A2, P) != 0 )
- return LW_FALSE;
-
- /* The point is on the segment line, check it's within its extent */
- return (P->x >= FP_MIN(A1->x, A2->x) && P->x <= FP_MAX(A1->x, A2->x) &&
- P->y >= FP_MIN(A1->y, A2->y) && P->y <= FP_MAX(A1->y, A2->y));
-}
-
-/**
-* How do the segment (P1,P2) and the segment (Q1,Q2) intersect?
-*
-* Only 2D coordinates are considered. Zero length segments are supported,
-* they can only intersect at their single point.
-*
-* @param P1 first point of the first segment
-* @param P2 second point of the first segment
-* @param Q1 first point of the second segment
-* @param Q2 second point of the second segment
-* @param ipoint storage for the intersection point, set when the returned
-* value is #LWSEG_INTERSECT_POINT
-*
-* @return one of the #LWSEGMENT_INTERSECTION values
-*/
-static int
-lwsegment_intersection(const POINT2D *P1,
- const POINT2D *P2,
- const POINT2D *Q1,
- const POINT2D *Q2,
- POINT2D *ipoint)
-{
- int s11, s12, s21, s22;
- int zeroP, zeroQ;
- double len2, tQ1, tQ2, tmin, tmax;
-
- zeroP = P2D_SAME_STRICT(P1, P2);
- zeroQ = P2D_SAME_STRICT(Q1, Q2);
-
- /* Zero length segments, if they intersect, do it on their single point */
- if ( zeroP || zeroQ )
- {
- if ( zeroP && zeroQ )
- {
- if ( ! P2D_SAME_STRICT(P1, Q1) )
- return LWSEG_INTERSECT_NONE;
- }
- else if ( zeroP )
- {
- if ( ! lw_pt_on_segment_inclusive(Q1, Q2, P1) )
- return LWSEG_INTERSECT_NONE;
- }
- else
- {
- if ( ! lw_pt_on_segment_inclusive(P1, P2, Q1) )
- return LWSEG_INTERSECT_NONE;
- }
-
- *ipoint = zeroP ? *P1 : *Q1;
- return LWSEG_INTERSECT_POINT;
- }
-
- /* Are both endpoints of the second segment on the same side of the first? */
- s11 = lw_segment_side(P1, P2, Q1);
- s12 = lw_segment_side(P1, P2, Q2);
- if ( (s11 > 0 && s12 > 0) || (s11 < 0 && s12 < 0) )
- return LWSEG_INTERSECT_NONE;
-
- /* Are both endpoints of the first segment on the same side of the second? */
- s21 = lw_segment_side(Q1, Q2, P1);
- s22 = lw_segment_side(Q1, Q2, P2);
- if ( (s21 > 0 && s22 > 0) || (s21 < 0 && s22 < 0) )
- return LWSEG_INTERSECT_NONE;
-
- /*
- * Both endpoints of the second segment are on the line of the first,
- * the segments are collinear and the intersection is their overlap.
- *
- * Use the projection on the (P1,P2) vector as the segment parameter,
- * the overlap is the [0,len2] interval trimmed by the [tQ1,tQ2] one.
- */
- if ( s11 == 0 && s12 == 0 )
- {
- len2 = (P2->x - P1->x) * (P2->x - P1->x) +
- (P2->y - P1->y) * (P2->y - P1->y);
- tQ1 = (Q1->x - P1->x) * (P2->x - P1->x) +
- (Q1->y - P1->y) * (P2->y - P1->y);
- tQ2 = (Q2->x - P1->x) * (P2->x - P1->x) +
- (Q2->y - P1->y) * (P2->y - P1->y);
-
- tmin = FP_MAX(0.0, FP_MIN(tQ1, tQ2));
- tmax = FP_MIN(len2, FP_MAX(tQ1, tQ2));
-
- /* No overlap at all */
- if ( tmin > tmax )
- return LWSEG_INTERSECT_NONE;
-
- /* The segments are collinear and share a positive length */
- if ( tmin < tmax )
- return LWSEG_INTERSECT_OVERLAP;
-
- /* The segments touch on a single point, always a segment endpoint */
- if ( tmin == 0.0 )
- *ipoint = *P1;
- else if ( tmin == len2 )
- *ipoint = *P2;
- else if ( tmin == tQ1 )
- *ipoint = *Q1;
- else if ( tmin == tQ2 )
- *ipoint = *Q2;
- else
- {
- /* Rounding fallback, the point is one of the four endpoints */
- double t = tmin / len2;
- ipoint->x = P1->x + t * (P2->x - P1->x);
- ipoint->y = P1->y + t * (P2->y - P1->y);
- }
- return LWSEG_INTERSECT_POINT;
- }
-
- /*
- * The segments are not collinear, they intersect on a single point.
- * When that point is not a vertex of either segment it is interior
- * to both of them.
- */
- if ( lw_pt_on_segment_inclusive(P1, P2, Q1) )
- {
- *ipoint = *Q1;
- return LWSEG_INTERSECT_POINT;
- }
-
- if ( lw_pt_on_segment_inclusive(P1, P2, Q2) )
- {
- *ipoint = *Q2;
- return LWSEG_INTERSECT_POINT;
- }
-
- if ( lw_pt_on_segment_inclusive(Q1, Q2, P1) )
- {
- *ipoint = *P1;
- return LWSEG_INTERSECT_POINT;
- }
-
- if ( lw_pt_on_segment_inclusive(Q1, Q2, P2) )
- {
- *ipoint = *P2;
- return LWSEG_INTERSECT_POINT;
- }
-
- return LWSEG_INTERSECT_CROSS;
-}
-
-/**
-* Is the point a boundary point of the line, that is one of its endpoints?
-*
-* Those are the only points of a line that may be shared with another line
-* without requiring any of them to be split.
-*
-* @param line the line to test the point against
-* @param P the point to test
-*
-* @return LW_TRUE if P is the first or the last point of the line
-*/
-static int
-lwline_point_is_boundary(const LWLINE *line, const POINT2D *P)
-{
- const POINT2D *start = getPoint2d_cp(line->points, 0);
- const POINT2D *end = getPoint2d_cp(line->points, line->points->npoints - 1);
-
- return P2D_SAME_STRICT(P, start) || P2D_SAME_STRICT(P, end);
-}
-
-/**
-* Do two lines intersect anywhere but on a point shared by both of their
-* boundaries?
-*
-* Only 2D coordinates are considered.
-*
-* @param line1 first line
-* @param line2 second line
-*
-* @return 1 if there's a non-boundary 2d intersection, 0 if there's NO
-* non-boundary 2d intersection
-*/
-int
-lwline_have_nonboundary_2d_intersection(const LWLINE *line1, const LWLINE *line2)
-{
- const POINTARRAY *pa1, *pa2;
- uint32_t i1, i2;
- POINT2D ipoint;
-
- if ( !line1 || !line2 )
- return LW_FALSE;
-
- pa1 = line1->points;
- pa2 = line2->points;
-
- /* Lines without segments can't intersect anything */
- if ( !pa1 || !pa2 || pa1->npoints < 2 || pa2->npoints < 2 )
- return LW_FALSE;
-
- for (i1 = 0; i1 < pa1->npoints - 1; ++i1)
- {
- const POINT2D *p1 = getPoint2d_cp(pa1, i1);
- const POINT2D *p2 = getPoint2d_cp(pa1, i1 + 1);
-
- for (i2 = 0; i2 < pa2->npoints - 1; ++i2)
- {
- const POINT2D *q1 = getPoint2d_cp(pa2, i2);
- const POINT2D *q2 = getPoint2d_cp(pa2, i2 + 1);
- int inter;
-
- inter = lwsegment_intersection(p1, p2, q1, q2, &ipoint);
-
- /* Crossings and overlaps are never boundary only */
- if ( inter == LWSEG_INTERSECT_CROSS ||
- inter == LWSEG_INTERSECT_OVERLAP )
- return LW_TRUE;
-
- /*
- ** A touch on a point which is a boundary point of both lines
- ** is not a non-boundary intersection
- */
- if ( inter == LWSEG_INTERSECT_POINT &&
- ( ! lwline_point_is_boundary(line1, &ipoint) ||
- ! lwline_point_is_boundary(line2, &ipoint) ) )
- return LW_TRUE;
- }
- }
-
- return LW_FALSE;
-}
diff --git a/liblwgeom/topo/liblwgeom_topo.h b/liblwgeom/topo/liblwgeom_topo.h
index fb25c3a10..24967e6aa 100644
--- a/liblwgeom/topo/liblwgeom_topo.h
+++ b/liblwgeom/topo/liblwgeom_topo.h
@@ -1432,4 +1432,26 @@ LWGEOM* lwt_GetFaceGeometry(LWT_TOPOLOGY* topo, LWT_ELEMID face);
int lwt_IsTopoRingCCW(const POINTARRAY *pa);
+/**
+ * Do two lines intersect anywhere but on a point shared by both of their
+ * boundaries, evaluated with exact arithmetic?
+ *
+ * Only 2D coordinates are considered. The classification of each pair of
+ * segments is evaluated with exact rational arithmetic, so that floating
+ * point cancellation between far-away coordinates cannot make a disjoint
+ * pair of lines look like an intersecting one. Single-point
+ * intersections are always a vertex of one of the two segments: the
+ * point is copied verbatim, never computed, so its comparison against
+ * the line boundary points is exact. Lines with non-finite
+ * coordinates (NaN or Inf) are reported as intersecting, so the caller
+ * rejects them as crossing.
+ *
+ * @param line1 first line
+ * @param line2 second line
+ *
+ * @return 1 if there's a non-boundary 2d intersection, 0 if there's NO
+ * non-boundary 2d intersection
+ */
+int lwt_LineHaveNonBoundary2DIntersection(const LWLINE *line1, const LWLINE *line2);
+
#endif /* LIBLWGEOM_TOPO_H */
diff --git a/liblwgeom/topo/lwgeom_topo.c b/liblwgeom/topo/lwgeom_topo.c
index 91efe0a20..3d8f7a8f3 100644
--- a/liblwgeom/topo/lwgeom_topo.c
+++ b/liblwgeom/topo/lwgeom_topo.c
@@ -625,6 +625,494 @@ lwt_AddIsoNode( LWT_TOPOLOGY* topo, LWT_ELEMID face,
return _lwt_AddIsoNode( topo, face, pt, skipISOChecks, 1 );
}
+/**
+* Exact sign of the 2D cross product (q-p1) x (p2-p1), that is the
+* orientation of q with respect to the directed segment (p1,p2).
+*
+* The computation is carried out with exact rational arithmetic using
+* the GNU Multiple Precision Library, so the sign cannot be fooled by
+* floating point cancellation when the magnitude of the coordinates is
+* much larger than the distances between the segments.
+*
+* The coordinates are expected to be finite, non-finite values are not
+* representable in the rationals.
+*
+* @return -1, 0 or 1, the sign of the exact cross product
+*/
+static int
+lwt_segment_side_exact(const POINT2D *p1, const POINT2D *p2, const POINT2D *q)
+{
+ mpq_t x1, y1, x2, y2, xq, yq;
+ mpq_t vx, vy, t1, t2, t3;
+ int side;
+
+ mpq_init(x1);
+ mpq_init(y1);
+ mpq_init(x2);
+ mpq_init(y2);
+ mpq_init(xq);
+ mpq_init(yq);
+ mpq_init(vx);
+ mpq_init(vy);
+ mpq_init(t1);
+ mpq_init(t2);
+ mpq_init(t3);
+
+ /* Exact conversion of the double coordinates to rationals */
+ mpq_set_d(x1, p1->x);
+ mpq_set_d(y1, p1->y);
+ mpq_set_d(x2, p2->x);
+ mpq_set_d(y2, p2->y);
+ mpq_set_d(xq, q->x);
+ mpq_set_d(yq, q->y);
+
+ /* v = p2 - p1 */
+ mpq_sub(vx, x2, x1);
+ mpq_sub(vy, y2, y1);
+
+ /* t3 = (qx-p1x) * (p2y-p1y) - (p2x-p1x) * (qy-p1y) */
+ mpq_sub(t1, xq, x1);
+ mpq_sub(t2, yq, y1);
+ mpq_mul(t3, t1, vy);
+ mpq_mul(t1, vx, t2);
+ mpq_sub(t3, t3, t1);
+
+ side = mpq_sgn(t3);
+
+ mpq_clear(x1);
+ mpq_clear(y1);
+ mpq_clear(x2);
+ mpq_clear(y2);
+ mpq_clear(xq);
+ mpq_clear(yq);
+ mpq_clear(vx);
+ mpq_clear(vy);
+ mpq_clear(t1);
+ mpq_clear(t2);
+ mpq_clear(t3);
+
+ return side;
+}
+
+/**
+* Is P in the extent (bounding box) of the (A1,A2) segment?
+*
+* All the comparisons are exact.
+*/
+static int
+lwt_pt_in_extent(const POINT2D *A1, const POINT2D *A2, const POINT2D *P)
+{
+ return (P->x >= FP_MIN(A1->x, A2->x) && P->x <= FP_MAX(A1->x, A2->x) &&
+ P->y >= FP_MIN(A1->y, A2->y) && P->y <= FP_MAX(A1->y, A2->y));
+}
+
+/**
+* Is P on the (A1,A2) segment, evaluated with exact arithmetic?
+*
+* A point is on the segment when its exact orientation with respect
+* to it is zero and it is within its extent.
+*/
+static int
+lwt_pt_on_segment_exact(const POINT2D *A1, const POINT2D *A2, const POINT2D *P)
+{
+ if ( lwt_segment_side_exact(A1, A2, P) != 0 )
+ return LW_FALSE;
+
+ /* The point is on the segment line, check it's within its extent */
+ return lwt_pt_in_extent(A1, A2, P);
+}
+
+/**
+* How do the segment (P1,P2) and the segment (Q1,Q2) intersect?
+*
+* All the orientation and projection tests are evaluated with exact
+* rational arithmetic, so the result is not affected by floating
+* point cancellation. Only finite coordinates are supported.
+*
+* The intersection point, when the segments intersect on a single
+* point, is never computed: any such point is a vertex of one of
+* the two segments, so it is copied verbatim from the segment
+* coordinates into ipoint. Copying an existing double does not
+* round anything, the point compares exactly against the vertex it
+* came from.
+*
+* @param P1 first point of the first segment
+* @param P2 second point of the first segment
+* @param Q1 first point of the second segment
+* @param Q2 second point of the second segment
+* @param ipoint storage for the intersection point, set when the returned
+* value is #LWSEG_INTERSECT_POINT
+*
+* @return one of the #LWSEGMENT_INTERSECTION values
+*/
+static int
+lwt_segment_intersection_exact(const POINT2D *P1,
+ const POINT2D *P2,
+ const POINT2D *Q1,
+ const POINT2D *Q2,
+ POINT2D *ipoint)
+{
+ mpq_t x1, y1, x2, y2, q1x, q1y, q2x, q2y;
+ mpq_t vx, vy, ux, uy;
+ mpq_t t1, t2, t3, t4;
+ mpq_t len2, tQ1, tQ2, tmin, tmax;
+ mpq_t zero;
+ int s11, s12, s21, s22;
+ int zeroP, zeroQ;
+ int ret;
+
+ zeroP = P2D_SAME_STRICT(P1, P2);
+ zeroQ = P2D_SAME_STRICT(Q1, Q2);
+
+ /* Zero length segments, if they intersect, do it on their single point */
+ if ( zeroP || zeroQ )
+ {
+ if ( zeroP && zeroQ )
+ {
+ if ( ! P2D_SAME_STRICT(P1, Q1) )
+ return LWSEG_INTERSECT_NONE;
+ }
+ else if ( zeroP )
+ {
+ if ( ! lwt_pt_on_segment_exact(Q1, Q2, P1) )
+ return LWSEG_INTERSECT_NONE;
+ }
+ else
+ {
+ if ( ! lwt_pt_on_segment_exact(P1, P2, Q1) )
+ return LWSEG_INTERSECT_NONE;
+ }
+
+ *ipoint = zeroP ? *P1 : *Q1;
+ return LWSEG_INTERSECT_POINT;
+ }
+
+ mpq_init(x1);
+ mpq_init(y1);
+ mpq_init(x2);
+ mpq_init(y2);
+ mpq_init(q1x);
+ mpq_init(q1y);
+ mpq_init(q2x);
+ mpq_init(q2y);
+ mpq_init(vx);
+ mpq_init(vy);
+ mpq_init(ux);
+ mpq_init(uy);
+ mpq_init(t1);
+ mpq_init(t2);
+ mpq_init(t3);
+ mpq_init(t4);
+ mpq_init(len2);
+ mpq_init(tQ1);
+ mpq_init(tQ2);
+ mpq_init(tmin);
+ mpq_init(tmax);
+ mpq_init(zero);
+ mpq_set_ui(zero, 0, 1);
+
+ /* Exact conversion of the double coordinates to rationals */
+ mpq_set_d(x1, P1->x);
+ mpq_set_d(y1, P1->y);
+ mpq_set_d(x2, P2->x);
+ mpq_set_d(y2, P2->y);
+ mpq_set_d(q1x, Q1->x);
+ mpq_set_d(q1y, Q1->y);
+ mpq_set_d(q2x, Q2->x);
+ mpq_set_d(q2y, Q2->y);
+
+ /* v = P2 - P1 */
+ mpq_sub(vx, x2, x1);
+ mpq_sub(vy, y2, y1);
+
+ /* s11 = sign( (Q1-P1) x (P2-P1) ) */
+ mpq_sub(t1, q1x, x1);
+ mpq_sub(t2, q1y, y1);
+ mpq_mul(t3, t1, vy);
+ mpq_mul(t4, vx, t2);
+ mpq_sub(t1, t3, t4);
+ s11 = mpq_sgn(t1);
+
+ /* s12 = sign( (Q2-P1) x (P2-P1) ) */
+ mpq_sub(t1, q2x, x1);
+ mpq_sub(t2, q2y, y1);
+ mpq_mul(t3, t1, vy);
+ mpq_mul(t4, vx, t2);
+ mpq_sub(t1, t3, t4);
+ s12 = mpq_sgn(t1);
+
+ /* Are both endpoints of the second segment on the same side of the first? */
+ if ( (s11 > 0 && s12 > 0) || (s11 < 0 && s12 < 0) )
+ {
+ ret = LWSEG_INTERSECT_NONE;
+ goto cleanup;
+ }
+
+ /* u = Q2 - Q1 */
+ mpq_sub(ux, q2x, q1x);
+ mpq_sub(uy, q2y, q1y);
+
+ /* s21 = sign( (P1-Q1) x (Q2-Q1) ) */
+ mpq_sub(t1, x1, q1x);
+ mpq_sub(t2, y1, q1y);
+ mpq_mul(t3, t1, uy);
+ mpq_mul(t4, ux, t2);
+ mpq_sub(t1, t3, t4);
+ s21 = mpq_sgn(t1);
+
+ /* s22 = sign( (P2-Q1) x (Q2-Q1) ) */
+ mpq_sub(t1, x2, q1x);
+ mpq_sub(t2, y2, q1y);
+ mpq_mul(t3, t1, uy);
+ mpq_mul(t4, ux, t2);
+ mpq_sub(t1, t3, t4);
+ s22 = mpq_sgn(t1);
+
+ /* Are both endpoints of the first segment on the same side of the second? */
+ if ( (s21 > 0 && s22 > 0) || (s21 < 0 && s22 < 0) )
+ {
+ ret = LWSEG_INTERSECT_NONE;
+ goto cleanup;
+ }
+
+ /*
+ * Both endpoints of the second segment are on the line of the first,
+ * the segments are collinear and the intersection is their overlap.
+ *
+ * Use the projection on the (P1,P2) vector as the segment parameter,
+ * the overlap is the [0,len2] interval trimmed by the [tQ1,tQ2] one.
+ * All the computations are exact.
+ */
+ if ( s11 == 0 && s12 == 0 )
+ {
+ mpq_mul(len2, vx, vx);
+ mpq_mul(t1, vy, vy);
+ mpq_add(len2, len2, t1);
+
+ /* tQ1 = (Q1-P1) . (P2-P1) */
+ mpq_sub(t1, q1x, x1);
+ mpq_mul(t2, t1, vx);
+ mpq_sub(t1, q1y, y1);
+ mpq_mul(t3, t1, vy);
+ mpq_add(tQ1, t2, t3);
+
+ /* tQ2 = (Q2-P1) . (P2-P1) */
+ mpq_sub(t1, q2x, x1);
+ mpq_mul(t2, t1, vx);
+ mpq_sub(t1, q2y, y1);
+ mpq_mul(t3, t1, vy);
+ mpq_add(tQ2, t2, t3);
+
+ /* tmin = max(0, min(tQ1, tQ2)) */
+ if ( mpq_cmp(tQ1, tQ2) <= 0 )
+ mpq_set(tmin, tQ1);
+ else
+ mpq_set(tmin, tQ2);
+ if ( mpq_cmp(tmin, zero) < 0 )
+ mpq_set(tmin, zero);
+
+ /* tmax = min(len2, max(tQ1, tQ2)) */
+ if ( mpq_cmp(tQ1, tQ2) > 0 )
+ mpq_set(tmax, tQ1);
+ else
+ mpq_set(tmax, tQ2);
+ if ( mpq_cmp(tmax, len2) > 0 )
+ mpq_set(tmax, len2);
+
+ /* No overlap at all */
+ if ( mpq_cmp(tmin, tmax) > 0 )
+ {
+ ret = LWSEG_INTERSECT_NONE;
+ goto cleanup;
+ }
+
+ /* The segments are collinear and share a positive length */
+ if ( mpq_cmp(tmin, tmax) < 0 )
+ {
+ ret = LWSEG_INTERSECT_OVERLAP;
+ goto cleanup;
+ }
+
+ /*
+ ** The segments touch on a single point. With non-zero-length
+ ** segments this is always one of the four segment endpoints:
+ ** tmin == tmax implies the touching parameter is 0, len2, tQ1
+ ** or tQ2, because otherwise tQ1 == tQ2 would make the second
+ ** segment zero-length. The point is copied verbatim from the
+ ** segment coordinates, never computed: no rounding of the
+ ** intersection point is involved.
+ */
+ if ( mpq_cmp(tmin, zero) == 0 )
+ *ipoint = *P1;
+ else if ( mpq_cmp(tmin, len2) == 0 )
+ *ipoint = *P2;
+ else if ( mpq_cmp(tmin, tQ1) == 0 )
+ *ipoint = *Q1;
+ else if ( mpq_cmp(tmin, tQ2) == 0 )
+ *ipoint = *Q2;
+ else
+ lwerror("%s: unreachable single-point collinear touch", __func__);
+ ret = LWSEG_INTERSECT_POINT;
+ goto cleanup;
+ }
+
+ /*
+ * The segments are not collinear, they intersect on a single point.
+ * When that point is not a vertex of either segment it is interior
+ * to both of them. The on-segment tests are exact.
+ */
+ if ( s11 == 0 && lwt_pt_in_extent(P1, P2, Q1) )
+ {
+ *ipoint = *Q1;
+ ret = LWSEG_INTERSECT_POINT;
+ goto cleanup;
+ }
+
+ if ( s12 == 0 && lwt_pt_in_extent(P1, P2, Q2) )
+ {
+ *ipoint = *Q2;
+ ret = LWSEG_INTERSECT_POINT;
+ goto cleanup;
+ }
+
+ if ( s21 == 0 && lwt_pt_in_extent(Q1, Q2, P1) )
+ {
+ *ipoint = *P1;
+ ret = LWSEG_INTERSECT_POINT;
+ goto cleanup;
+ }
+
+ if ( s22 == 0 && lwt_pt_in_extent(Q1, Q2, P2) )
+ {
+ *ipoint = *P2;
+ ret = LWSEG_INTERSECT_POINT;
+ goto cleanup;
+ }
+
+ ret = LWSEG_INTERSECT_CROSS;
+
+cleanup:
+ mpq_clear(x1);
+ mpq_clear(y1);
+ mpq_clear(x2);
+ mpq_clear(y2);
+ mpq_clear(q1x);
+ mpq_clear(q1y);
+ mpq_clear(q2x);
+ mpq_clear(q2y);
+ mpq_clear(vx);
+ mpq_clear(vy);
+ mpq_clear(ux);
+ mpq_clear(uy);
+ mpq_clear(t1);
+ mpq_clear(t2);
+ mpq_clear(t3);
+ mpq_clear(t4);
+ mpq_clear(len2);
+ mpq_clear(tQ1);
+ mpq_clear(tQ2);
+ mpq_clear(tmin);
+ mpq_clear(tmax);
+ mpq_clear(zero);
+
+ return ret;
+}
+
+/*
+ * Do two lines intersect anywhere but on a point shared by both of their
+ * boundaries, evaluated with exact arithmetic?
+ *
+ * Only 2D coordinates are considered.
+ *
+ * All the orientation and projection tests are evaluated with exact
+ * rational arithmetic (see lwt_segment_intersection_exact above), so that
+ * floating point cancellation between far-away coordinates cannot make a
+ * disjoint pair of lines look like an intersecting one.
+ *
+ * Coordinates that are not finite (NaN or Inf) cannot be converted to
+ * rationals; lines with such coordinates are reported as having a
+ * non-boundary intersection, which makes the caller reject them as
+ * crossing.
+ *
+ * @param line1 first line
+ * @param line2 second line
+ *
+ * @return 1 if there's a non-boundary 2d intersection, 0 if there's NO
+ * non-boundary 2d intersection
+ */
+int
+lwt_LineHaveNonBoundary2DIntersection(const LWLINE *line1, const LWLINE *line2)
+{
+ const POINTARRAY *pa1, *pa2;
+ uint32_t i1, i2;
+ POINT2D ipoint;
+ const POINT2D *s1, *e1, *s2, *e2;
+
+ if ( !line1 || !line2 )
+ return LW_FALSE;
+
+ pa1 = line1->points;
+ pa2 = line2->points;
+
+ /* Lines without segments can't intersect anything */
+ if ( !pa1 || !pa2 || pa1->npoints < 2 || pa2->npoints < 2 )
+ return LW_FALSE;
+
+ /* Non-finite coordinates cannot be converted to rationals: report
+ * them as intersecting, so the caller rejects the line as crossing */
+ for (i1 = 0; i1 < pa1->npoints; ++i1)
+ {
+ const POINT2D *p = getPoint2d_cp(pa1, i1);
+ if ( ! isfinite(p->x) || ! isfinite(p->y) )
+ return LW_TRUE;
+ }
+ for (i2 = 0; i2 < pa2->npoints; ++i2)
+ {
+ const POINT2D *p = getPoint2d_cp(pa2, i2);
+ if ( ! isfinite(p->x) || ! isfinite(p->y) )
+ return LW_TRUE;
+ }
+
+ /* Boundary points of the two lines */
+ s1 = getPoint2d_cp(pa1, 0);
+ e1 = getPoint2d_cp(pa1, pa1->npoints - 1);
+ s2 = getPoint2d_cp(pa2, 0);
+ e2 = getPoint2d_cp(pa2, pa2->npoints - 1);
+
+ for (i1 = 0; i1 < pa1->npoints - 1; ++i1)
+ {
+ const POINT2D *p1 = getPoint2d_cp(pa1, i1);
+ const POINT2D *p2 = getPoint2d_cp(pa1, i1 + 1);
+
+ for (i2 = 0; i2 < pa2->npoints - 1; ++i2)
+ {
+ const POINT2D *q1 = getPoint2d_cp(pa2, i2);
+ const POINT2D *q2 = getPoint2d_cp(pa2, i2 + 1);
+ int inter;
+
+ inter = lwt_segment_intersection_exact(p1, p2, q1, q2, &ipoint);
+
+ /* Crossings and overlaps are never boundary only */
+ if ( inter == LWSEG_INTERSECT_CROSS ||
+ inter == LWSEG_INTERSECT_OVERLAP )
+ return LW_TRUE;
+
+ /*
+ ** A touch on a point which is a boundary point of both lines
+ ** is not a non-boundary intersection
+ */
+ if ( inter == LWSEG_INTERSECT_POINT &&
+ ( ! (P2D_SAME_STRICT(&ipoint, s1) || P2D_SAME_STRICT(&ipoint, e1)) ||
+ ! (P2D_SAME_STRICT(&ipoint, s2) || P2D_SAME_STRICT(&ipoint, e2)) ) )
+ return LW_TRUE;
+ }
+ }
+
+ return LW_FALSE;
+}
+
+
+
/*
* Check that an edge does not cross an existing node and
* does not have non-boundary intersection with existing edge
@@ -710,7 +1198,7 @@ _lwt_CheckEdgeCrossing( LWT_TOPOLOGY* topo,
}
/* check if the edge has a non-boundary-boundary intersection with our edge */
- if ( lwline_have_nonboundary_2d_intersection(edge->geom, geom) )
+ if ( lwt_LineHaveNonBoundary2DIntersection(edge->geom, geom) )
{
char locinfo[128] = "";
_lwt_describe_intersection_point(lwline_as_lwgeom(edge->geom), lwline_as_lwgeom(geom), locinfo, sizeof(locinfo));
diff --git a/topology/postgis_topology.c b/topology/postgis_topology.c
index c96637b6e..43b5ad25e 100644
--- a/topology/postgis_topology.c
+++ b/topology/postgis_topology.c
@@ -5988,7 +5988,7 @@ Datum TopoEdgesHaveInteriorIntersection(PG_FUNCTION_ARGS)
PG_RETURN_NULL();
}
- ret = lwline_have_nonboundary_2d_intersection(lwgeom_as_lwline(lwgeom1), lwgeom_as_lwline(lwgeom2));
+ ret = lwt_LineHaveNonBoundary2DIntersection(lwgeom_as_lwline(lwgeom1), lwgeom_as_lwline(lwgeom2));
lwgeom_free(lwgeom1);
lwgeom_free(lwgeom2);
diff --git a/topology/test/regress/st_addedgemodface.sql b/topology/test/regress/st_addedgemodface.sql
index fdbb4dcf4..6384adb4c 100644
--- a/topology/test/regress/st_addedgemodface.sql
+++ b/topology/test/regress/st_addedgemodface.sql
@@ -556,6 +556,40 @@ SELECT NULL FROM ST_AddEdgeModFace('t6140', 1, 3, '0102000000030000005163B458843
SELECT 't6140', 'unexpected success', * FROM ST_AddEdgeModFace('t6140', 2, 3, '010200000002000000685E205B88333540B7363CF4E66A5140831D215B88333540107B41F4E66A5140');
ROLLBACK;
+-----------------------------------------------------
+-- Very far-away but disjoint edges
+-----------------------------------------------------
+
+-- See https://trac.osgeo.org/postgis/ticket/6143
+-- The two edges are disjoint: with exact arithmetic the distance
+-- between them is about 0.000070710678, so GEOS computes an empty
+-- intersection for the pair. The floating point implementation of
+-- the internal edge crossing check used to report an interior
+-- intersection, making ST_AddEdgeModFace reject the (valid) second
+-- edge with
+--
+-- SQL/MM Spatial exception - non-boundary intersection
+-- with edge 1 at or near
+--
+-- and, on commit d03570c85 (first implementation of the internal
+-- crossing check), crashing the backend while computing that
+-- "at or near" location.
+--
+-- The predicate now evaluates the intersections with exact rational
+-- arithmetic, so the correct behaviour (edge 2 accepted) is what
+-- this test expects.
+BEGIN;
+SELECT 't6143', 'start';
+SELECT NULL FROM CreateTopology('t6143');
+SELECT NULL FROM ST_AddIsoNode('t6143', NULL, 'POINT(-1000000000000 -1000000000000)'); -- start node of edge 1
+SELECT NULL FROM ST_AddIsoNode('t6143', NULL, 'POINT(-1 0.0001)'); -- end node of edge 1
+SELECT NULL FROM ST_AddIsoEdge('t6143', 1, 2, 'LINESTRING(-1000000000000 -1000000000000, -1 0.0001)');
+SELECT NULL FROM ST_AddIsoNode('t6143', NULL, 'POINT(100000 -10000)'); -- start node of edge 2
+SELECT NULL FROM ST_AddIsoNode('t6143', NULL, 'POINT(-1 -0.00000000001)'); -- end node of edge 2
+-- This is the edge that used to be wrongly rejected
+SELECT 't6143', 'edge accepted', ST_AddEdgeModFace('t6143', 3, 4, 'LINESTRING(100000 -10000, -1 -0.00000000001)');
+ROLLBACK;
+
---------------------------------------------------------------------
-- Cleanups
---------------------------------------------------------------------
diff --git a/topology/test/regress/st_addedgemodface_expected b/topology/test/regress/st_addedgemodface_expected
index 4b92f3a5f..49a0bc322 100644
--- a/topology/test/regress/st_addedgemodface_expected
+++ b/topology/test/regress/st_addedgemodface_expected
@@ -205,4 +205,6 @@ F35|POLYGON((9 14,9 22,21 22,21 14,9 14))
ERROR: Backend error: Corrupted topology: ring of edge 58 is topologically non-closed
t6140|start
ERROR: SQL/MM Spatial exception - non-boundary intersection with edge 1 at or near POINT(21.2012993770454 69.6703463210545)
+t6143|start
+t6143|edge accepted|2
Topology 'city_data' dropped
-----------------------------------------------------------------------
Summary of changes:
liblwgeom/cunit/cu_algorithm.c | 153 -------
liblwgeom/cunit/cu_topo.c | 161 ++++++++
liblwgeom/liblwgeom_internal.h | 18 +-
liblwgeom/lwline.c | 261 ------------
liblwgeom/topo/liblwgeom_topo.h | 22 +
liblwgeom/topo/lwgeom_topo.c | 490 ++++++++++++++++++++++-
topology/postgis_topology.c | 2 +-
topology/test/regress/st_addedgemodface.sql | 34 ++
topology/test/regress/st_addedgemodface_expected | 2 +
9 files changed, 726 insertions(+), 417 deletions(-)
hooks/post-receive
--
PostGIS
More information about the postgis-tickets
mailing list