Fix intersection computation of slightly overlapping polygons

src/edgegraph/HalfEdge.cpp

@@ -179,17 +179,14 @@ HalfEdge::findLowest() const
 int
 HalfEdge::compareAngularDirection(const HalfEdge* e) const
 {
-    double dx = directionX();
-    double dy = directionY();
-    double dx2 = e->directionX();
-    double dy2 = e->directionY();
+    assert(orig().equals2D(e->orig()));
 
     // same vector
-    if (dx == dx2 && dy == dy2)
+    if (directionPt().equals2D(e->directionPt()))
         return 0;
 
-    int quadrant = geom::Quadrant::quadrant(dx, dy);
-    int quadrant2 = geom::Quadrant::quadrant(dx2, dy2);
+    int quadrant = geom::Quadrant::quadrant(directionX(), directionY());
+    int quadrant2 = geom::Quadrant::quadrant(e->directionX(), e->directionY());
 
     /**
     * If the direction vectors are in different quadrants,
@@ -272,5 +269,3 @@ HalfEdge::toStringNode(const HalfEdge* he, std::ostream& os)
 
 } // namespace geos.edgegraph
 } // namespace geos
-
-

src/geomgraph/EdgeEnd.cpp

@@ -160,7 +160,9 @@ int
 EdgeEnd::compareDirection(const EdgeEnd* e) const
 {
     assert(e);
-    if(dx == e->dx && dy == e->dy) {
+    assert(p0.equals2D(e->p0));
+
+    if(p1.equals2D(e->p1)) {
         return 0;
     }
 
@@ -215,4 +217,3 @@ operator<< (std::ostream& os, const EdgeEnd& ee)
 
 } // namespace geos.geomgraph
 } // namespace geos
-

tests/unit/edgegraph/EdgeGraphTest.cpp

@@ -170,6 +170,25 @@ void object::test<5> ()
     checkNextPrev(*graph);
 }
 
+// testSimilarEdgesDirection
+template<>
+template<>
+void object::test<6> ()
+{
+    EdgeGraph graph;
+    HalfEdge* e1 = addEdge(graph, 1, 1, 0, 0.5);
+    addEdge(graph, 1, 1, 0, 0.49999999999999994);
+    addEdge(graph, 1, 1, 0, 1);
+
+    HalfEdge* eUpper = findEdge(graph, 1, 1, 0, 0.5);
+    HalfEdge* eLower = findEdge(graph, 1, 1, 0, 0.49999999999999994);
+    ensure("edges with distinct direction points must not compare equal", eUpper->compareTo(eLower) != 0);
+    ensure("edge comparison must be antisymmetric", eUpper->compareTo(eLower) == -eLower->compareTo(eUpper));
+
+    checkNodeValid(e1);
+    checkNextPrev(graph);
+}
+
 
 
 } // namespace tut

tests/unit/operation/overlayng/OverlayNGRobustTest.cpp

@@ -61,6 +61,17 @@ struct test_overlayngrobust_data {
         ensure_NO_THROW( OverlayNGRobust::Overlay(geom_a.get(), geom_b.get(), opCode) );
     }
 
+    void
+    checkOverlaySymmetricArea(const std::string& a, const std::string& b, int opCode, double expectedArea)
+    {
+        std::unique_ptr<Geometry> geom_a = r.read(a);
+        std::unique_ptr<Geometry> geom_b = r.read(b);
+        double areaAB = OverlayNGRobust::Overlay(geom_a.get(), geom_b.get(), opCode)->getArea();
+        double areaBA = OverlayNGRobust::Overlay(geom_b.get(), geom_a.get(), opCode)->getArea();
+        ensure_equals("overlay is not commutative", areaAB, areaBA, 1e-12);
+        ensure_equals("unexpected overlay area", areaAB, expectedArea, 1e-12);
+    }
+
     std::unique_ptr<Geometry>
     double2geom(const std::vector<double>& x, const std::vector<double>& y)
     {
@@ -109,7 +120,26 @@ void object::test<2> ()
     checkOverlaySuccess(a, b, OverlayNG::INTERSECTION);
 }
 
+// 2026-08-05 Intersection of polygons isn't commutative (https://github.com/libgeos/geos/issues/1405)
+template<>
+template<>
+void object::test<3> ()
+{
+    set_test_name("Intersection of triangles with an extremely similar edge");
+    const std::string a = "POLYGON ((1 1, 0 0.5, 0 0, 1 1))";
+    const std::string b = "POLYGON ((1 1, 0 0.49999999999999994, 0 1, 1 1))";
+    checkOverlaySymmetricArea(a, b, OverlayNG::INTERSECTION, 0.0);
+}
 
+template<>
+template<>
+void object::test<4> ()
+{
+    set_test_name("Union of triangles with an extremely similar edge");
+    const std::string a = "POLYGON ((1 1, 0 0.5, 0 0, 1 1))";
+    const std::string b = "POLYGON ((1 1, 0 0.49999999999999994, 0 1, 1 1))";
+    checkOverlaySymmetricArea(a, b, OverlayNG::UNION, 0.5);
+}
 
 #if 0
 /**