From 91857fdabbd383f26635a4b9ac103777bee0a872 Mon Sep 17 00:00:00 2001
From: tanelk <tanel.kiis@gmail.com>
Date: Sun, 17 Mar 2024 17:16:34 +0200
Subject: [PATCH 01/12] Fix vertex-vertex intersection in OverlayArea

Signed-off-by: tanelk <tanel.kiis@gmail.com>
---
 .../operation/overlayarea/OverlayArea.java    | 204 ++++++++++--------
 .../overlayarea/BaseOverlayAreaTest.java      |  84 ++++++++
 .../overlayarea/OverlayAreaTest.java          |  60 +-----
 .../overlayarea/SimpleOverlayAreaTest.java    |  63 +-----
 4 files changed, 211 insertions(+), 200 deletions(-)
 create mode 100644 modules/lab/src/test/java/org/locationtech/jts/operation/overlayarea/BaseOverlayAreaTest.java

diff --git a/modules/lab/src/main/java/org/locationtech/jts/operation/overlayarea/OverlayArea.java b/modules/lab/src/main/java/org/locationtech/jts/operation/overlayarea/OverlayArea.java
index 8e475ebd52..58758ead8f 100644
--- a/modules/lab/src/main/java/org/locationtech/jts/operation/overlayarea/OverlayArea.java
+++ b/modules/lab/src/main/java/org/locationtech/jts/operation/overlayarea/OverlayArea.java
@@ -11,8 +11,7 @@
  */
 package org.locationtech.jts.operation.overlayarea;
 
-import java.util.List;
-
+import org.locationtech.jts.algorithm.Angle;
 import org.locationtech.jts.algorithm.Area;
 import org.locationtech.jts.algorithm.LineIntersector;
 import org.locationtech.jts.algorithm.Orientation;
@@ -24,16 +23,23 @@
 import org.locationtech.jts.geom.CoordinateSequence;
 import org.locationtech.jts.geom.Envelope;
 import org.locationtech.jts.geom.Geometry;
+import org.locationtech.jts.geom.GeometryFactory;
 import org.locationtech.jts.geom.GeometryFilter;
-import org.locationtech.jts.geom.LineSegment;
 import org.locationtech.jts.geom.LinearRing;
 import org.locationtech.jts.geom.Location;
+import org.locationtech.jts.geom.Point;
 import org.locationtech.jts.geom.Polygon;
-import org.locationtech.jts.index.ItemVisitor;
 import org.locationtech.jts.index.kdtree.KdNode;
 import org.locationtech.jts.index.kdtree.KdTree;
-import org.locationtech.jts.index.strtree.STRtree;
 import org.locationtech.jts.math.MathUtil;
+import org.locationtech.jts.noding.BasicSegmentString;
+import org.locationtech.jts.noding.MCIndexSegmentSetMutualIntersector;
+import org.locationtech.jts.noding.SegmentIntersector;
+import org.locationtech.jts.noding.SegmentSetMutualIntersector;
+import org.locationtech.jts.noding.SegmentString;
+
+import java.util.Collections;
+import java.util.List;
 
 /**
  * Computes the area of the overlay of two polygons without forming
@@ -79,7 +85,7 @@ private static boolean interacts(Geometry geom0, Geometry geom1) {
   private Geometry geom0;
   private Envelope geomEnv0;
   private IndexedPointInAreaLocator locator0;
-  private STRtree indexSegs;
+  private SegmentSetMutualIntersector segSetMutInt;
   private KdTree vertexIndex;
 
   public OverlayArea(Geometry geom) {
@@ -92,7 +98,7 @@ public OverlayArea(Geometry geom) {
     
     geomEnv0 = geom.getEnvelopeInternal();
     locator0 = new IndexedPointInAreaLocator(geom);
-    indexSegs = buildSegmentIndex(geom);
+    segSetMutInt = buildSegmentIndex(geom);
     vertexIndex = buildVertexIndex(geom);
   }
   
@@ -207,80 +213,114 @@ private static double area(Geometry geom) {
   }
 
   private double areaForIntersections(LinearRing geom) {
-    double area = 0.0;
-    CoordinateSequence seq = geom.getCoordinateSequence();
-    
-    boolean isCCW = Orientation.isCCW(seq);
-    
-    // Compute rays for all intersections   
-    for (int j = 0; j < seq.size() - 1; j++) {
-      Coordinate b0 = seq.getCoordinate(j);
-      Coordinate b1 = seq.getCoordinate(j+1);
-      if (isCCW) {
-        // flip segment orientation
-        Coordinate temp = b0; b0 = b1; b1 = temp;
-      }
-      
-      Envelope env = new Envelope(b0, b1);
-      IntersectionVisitor intVisitor = new IntersectionVisitor(b0, b1);
-      indexSegs.query(env, intVisitor);
-      area += intVisitor.getArea();
-    }
-    return area;
+    Coordinate[] coords = geom.getCoordinates();
+    SegmentString segStr = new BasicSegmentString(coords, Orientation.isCCW(coords));
+
+    IntersectionVisitor intVisitor = new IntersectionVisitor();
+    segSetMutInt.process(Collections.singletonList(segStr), intVisitor);
+
+    return intVisitor.getArea();
   }
 
-  class IntersectionVisitor implements ItemVisitor {
-    double area = 0.0;
-    private Coordinate b0;
-    private Coordinate b1;
-    
-    IntersectionVisitor(Coordinate b0, Coordinate b1) {
-      this.b0 = b0;
-      this.b1 = b1;
-    }
-    
+  class IntersectionVisitor implements SegmentIntersector {
+    private double area = 0.0;
+
     double getArea() {
       return area;
     }
-    
-    public void visitItem(Object item) {
-      LineSegment seg = (LineSegment) item;
-      area += areaForIntersection(seg.p0, seg.p1, b0, b1);
-    }
-  }
-  
-  private static double areaForIntersection(Coordinate a0, Coordinate a1, Coordinate b0, Coordinate b1 ) {
-    // TODO: can the intersection computation be optimized?
-    li.computeIntersection(a0, a1, b0, b1);
-    if (! li.hasIntersection()) return 0.0;
-    
-    /**
-     * An intersection creates two edge vectors which contribute to the area.
-     * 
-     * With both rings oriented CW (effectively)
-     * There are two situations for segment intersection:
-     * 
-     * 1) A entering B, B exiting A => rays are IP->A1:R, IP->B0:L
-     * 2) A exiting B, B entering A => rays are IP->A0:L, IP->B1:R
-     * (where IP is the intersection point, 
-     * and  :L/R indicates result polygon interior is to the Left or Right).
-     * 
-     * For accuracy the full edge is used to provide the direction vector.
-     */
-    Coordinate intPt = li.getIntersection(0);
-    
-    boolean isAenteringB = Orientation.COUNTERCLOCKWISE == Orientation.index(a0, a1, b1);
-    
-    if ( isAenteringB ) {
-      return EdgeVector.area2Term(intPt, a0, a1, true)
-        + EdgeVector.area2Term(intPt, b1, b0, false);
+
+    @Override
+    public void processIntersections(SegmentString a, int aIndex, SegmentString b, int bIndex) {
+      boolean isCCWA = (boolean) a.getData();
+      boolean isCCWB = (boolean) b.getData();
+
+      Coordinate a0 = a.getCoordinate(aIndex);
+      Coordinate a1 = a.getCoordinate(aIndex + 1);
+      Coordinate b0 = b.getCoordinate(bIndex);
+      Coordinate b1 = b.getCoordinate(bIndex + 1);
+
+      if (isCCWA) {
+        Coordinate tmp = a0; a0 = a1; a1 = tmp;
+      }
+      if (isCCWB) {
+          Coordinate tmp = b0; b0 = b1; b1 = tmp;
+      }
+
+      li.computeIntersection(a0, a1, b0, b1);
+      if (! li.hasIntersection()) return;
+
+      Coordinate intPt = li.getIntersection(0);
+
+      if (li.isProper() || li.isInteriorIntersection()) {
+        // Edge-edge intersection OR vertex-edge intersection
+
+        /**
+         * An intersection creates two edge vectors which contribute to the area.
+         *
+         * With both rings oriented CW (effectively)
+         * There are two situations for segment intersection:
+         *
+         * 1) A entering B, B exiting A => rays are IP->A1:R, IP->B0:L
+         * 2) A exiting B, B entering A => rays are IP->A0:L, IP->B1:R
+         * (where IP is the intersection point,
+         * and  :L/R indicates result polygon interior is to the Left or Right).
+         *
+         * For accuracy the full edge is used to provide the direction vector.
+         */
+
+        if (Orientation.CLOCKWISE == Orientation.index(a0, a1, b0)) {
+          area += EdgeVector.area2Term(intPt, a0, a1, true);
+          area += EdgeVector.area2Term(intPt, b1, b0, false);
+        } else if (Orientation.CLOCKWISE == Orientation.index(a0, a1, b1)) {
+          area += EdgeVector.area2Term(intPt, a1, a0, false);
+          area += EdgeVector.area2Term(intPt, b0, b1, true);
+        }
+
+      } else {
+        // vertex-vertex intersection
+        // This intersection is visited 4 times - include only once
+        if (!intPt.equals2D(a1) || !intPt.equals2D(b1)) {
+          return;
+        }
+
+        Coordinate a2 = a.nextInRing(aIndex + 1);
+        Coordinate b2 = b.nextInRing(bIndex + 1);
+        if (isCCWA) {
+          a2 = a.prevInRing(aIndex);
+        }
+        if (isCCWB) {
+          b2 = b.prevInRing(bIndex);
+        }
+
+        double aAngle = Angle.interiorAngle(a0, intPt, a2);
+        double bAngle = Angle.interiorAngle(b0, intPt, b2);
+
+        // The LTE ja LT are chosen such that when A0->A1 is collinear with B0->B1, then A is chosen
+        // and when A1->A2 is collinear with B1->B2, then A is chosen.
+        // This avoids double counting in the case of collinear segments.
+        if (Angle.interiorAngle(a0, intPt, b2) <= bAngle) {
+          area += EdgeVector.area2Term(intPt, a1, a0, false);
+        }
+        if (Angle.interiorAngle(b0, intPt, a2) <= bAngle) {
+          area += EdgeVector.area2Term(intPt, a1, a2, true);
+        }
+        if (Angle.interiorAngle(b0, intPt, a2) < aAngle) {
+          area += EdgeVector.area2Term(intPt, b1, b0, false);
+        }
+        if (Angle.interiorAngle(a0, intPt, b2) < aAngle) {
+          area += EdgeVector.area2Term(intPt, b1, b2, true);
+        }
+      }
     }
-    else {
-      return EdgeVector.area2Term(intPt, a1, a0, false)
-       + EdgeVector.area2Term(intPt, b0, b1, true);
+
+    @Override
+    public boolean isDone() {
+      // Process all intersections
+      return false;
     }
   }
-    
+  
+
   private double areaForInteriorVertices(LinearRing ring) {
     /**
      * Compute rays originating at vertices inside the intersection result
@@ -335,22 +375,10 @@ private static CoordinateSequence getVertices(Geometry geom) {
     return seq;
   }
   
-  private static STRtree buildSegmentIndex(Geometry geom) {
+  private static SegmentSetMutualIntersector buildSegmentIndex(Geometry geom) {
     Coordinate[] coords = geom.getCoordinates();
-    
-    boolean isCCW = Orientation.isCCW(coords);
-    STRtree index = new STRtree();
-    for (int i = 0; i < coords.length - 1; i++) {
-      Coordinate a0 = coords[i];
-      Coordinate a1 = coords[i+1];
-      LineSegment seg = new LineSegment(a0, a1);
-      if (isCCW) {
-        seg = new LineSegment(a1, a0);
-      }
-      Envelope env = new Envelope(a0, a1);
-      index.insert(env, seg);
-    }
-    return index;
+    SegmentString segStr = new BasicSegmentString(coords, Orientation.isCCW(coords));
+    return new MCIndexSegmentSetMutualIntersector(Collections.singletonList(segStr));
   }
 
   private static KdTree buildVertexIndex(Geometry geom) {
diff --git a/modules/lab/src/test/java/org/locationtech/jts/operation/overlayarea/BaseOverlayAreaTest.java b/modules/lab/src/test/java/org/locationtech/jts/operation/overlayarea/BaseOverlayAreaTest.java
new file mode 100644
index 0000000000..01feff1ea0
--- /dev/null
+++ b/modules/lab/src/test/java/org/locationtech/jts/operation/overlayarea/BaseOverlayAreaTest.java
@@ -0,0 +1,84 @@
+/*
+ * Copyright (c) 2022 Martin Davis
+ *
+ * All rights reserved. This program and the accompanying materials
+ * are made available under the terms of the Eclipse Public License 2.0
+ * and Eclipse Distribution License v. 1.0 which accompanies this distribution.
+ * The Eclipse Public License is available at http://www.eclipse.org/legal/epl-v20.html
+ * and the Eclipse Distribution License is available at
+ *
+ * http://www.eclipse.org/org/documents/edl-v10.php.
+ */
+package org.locationtech.jts.operation.overlayarea;
+
+import org.locationtech.jts.geom.Geometry;
+import test.jts.GeometryTestCase;
+
+/**
+ * Base class for testing overlay area algorithms
+ * by comparing them to the standard overlay intersection area.
+ * This is a parameterized test class.
+ * Subclasses provide the overlay area algorithm to test.
+ */
+public abstract class BaseOverlayAreaTest extends GeometryTestCase {
+
+  public BaseOverlayAreaTest(String name) {
+    super(name);
+  }
+
+  public void testDisjoint() {
+    checkIntersectionArea(
+        "POLYGON ((10 90, 40 90, 40 60, 10 60, 10 90))",
+        "POLYGON ((90 10, 50 10, 50 50, 90 50, 90 10))");
+  }
+
+  public void testTouching() {
+    checkIntersectionArea(
+        "POLYGON ((10 90, 50 90, 50 50, 10 50, 10 90))",
+        "POLYGON ((90 10, 50 10, 50 50, 90 50, 90 10))");
+  }
+
+  public void testRectangleAContainsB() {
+    checkIntersectionArea(
+        "POLYGON ((100 300, 300 300, 300 100, 100 100, 100 300))",
+        "POLYGON ((150 250, 250 250, 250 150, 150 150, 150 250))");
+  }
+
+  public void testTriangleAContainsB() {
+    checkIntersectionArea(
+        "POLYGON ((60 170, 270 370, 380 60, 60 170))",
+        "POLYGON ((200 250, 245 155, 291 195, 200 250))");
+  }
+
+  public void testRectangleOverlap() {
+    checkIntersectionArea(
+        "POLYGON ((100 200, 200 200, 200 100, 100 100, 100 200))",
+        "POLYGON ((250 250, 250 150, 150 150, 150 250, 250 250))");
+  }
+
+  public void testRectangleTriangleOverlap() {
+    checkIntersectionArea(
+        "POLYGON ((100 200, 200 200, 200 100, 100 100, 100 200))",
+        "POLYGON ((300 200, 150 150, 300 100, 300 200))");
+  }
+
+  public void testSawOverlap() {
+    checkIntersectionArea(
+        "POLYGON ((100 300, 305 299, 150 200, 300 150, 150 100, 300 50, 100 50, 100 300))",
+        "POLYGON ((400 350, 150 250, 350 200, 200 150, 350 100, 180 50, 400 50, 400 350))");
+  }
+
+  protected void checkIntersectionArea(String wktA, String wktB) {
+    Geometry a = read(wktA);
+    Geometry b = read(wktB);
+    
+    double ovIntArea = computeOverlayArea(a, b);
+    
+    double intAreaFull = a.intersection(b).getArea();
+    
+    //System.out.printf("OverlayArea: %f   Full overlay: %f\n", ovIntArea, intAreaFull);
+    assertEquals(intAreaFull, ovIntArea, 0.0001);
+  }
+
+  abstract protected double computeOverlayArea(Geometry a, Geometry b);
+}
diff --git a/modules/lab/src/test/java/org/locationtech/jts/operation/overlayarea/OverlayAreaTest.java b/modules/lab/src/test/java/org/locationtech/jts/operation/overlayarea/OverlayAreaTest.java
index 369ee013aa..8d79c9a17b 100644
--- a/modules/lab/src/test/java/org/locationtech/jts/operation/overlayarea/OverlayAreaTest.java
+++ b/modules/lab/src/test/java/org/locationtech/jts/operation/overlayarea/OverlayAreaTest.java
@@ -14,9 +14,10 @@
 import org.locationtech.jts.geom.Geometry;
 
 import junit.textui.TestRunner;
+import org.locationtech.jts.geom.Polygon;
 import test.jts.GeometryTestCase;
 
-public class OverlayAreaTest extends GeometryTestCase {
+public class OverlayAreaTest extends BaseOverlayAreaTest {
 
   public static void main(String args[]) {
     TestRunner.run(OverlayAreaTest.class);
@@ -26,49 +27,6 @@ public OverlayAreaTest(String name) {
     super(name);
   }
 
-  public void testDisjoint() {
-    checkIntersectionArea(
-        "POLYGON ((10 90, 40 90, 40 60, 10 60, 10 90))",
-        "POLYGON ((90 10, 50 10, 50 50, 90 50, 90 10))");
-  }
-  
-  //TODO: fix this bug
-  public void xtestTouching() {
-    checkIntersectionArea(
-        "POLYGON ((10 90, 50 90, 50 50, 10 50, 10 90))",
-        "POLYGON ((90 10, 50 10, 50 50, 90 50, 90 10))");
-  }
-  
-  public void testRectangleAContainsB() {
-    checkIntersectionArea(
-        "POLYGON ((100 300, 300 300, 300 100, 100 100, 100 300))",
-        "POLYGON ((150 250, 250 250, 250 150, 150 150, 150 250))");
-  }
-
-  public void testTriangleAContainsB() {
-    checkIntersectionArea(
-        "POLYGON ((60 170, 270 370, 380 60, 60 170))",
-        "POLYGON ((200 250, 245 155, 291 195, 200 250))");
-  }
-
-  public void testRectangleOverlap() {
-    checkIntersectionArea(
-        "POLYGON ((100 200, 200 200, 200 100, 100 100, 100 200))",
-        "POLYGON ((250 250, 250 150, 150 150, 150 250, 250 250))");
-  }
-
-  public void testRectangleTriangleOverlap() {
-    checkIntersectionArea(
-        "POLYGON ((100 200, 200 200, 200 100, 100 100, 100 200))",
-        "POLYGON ((300 200, 150 150, 300 100, 300 200))");
-  }
-
-  public void testSawOverlap() {
-    checkIntersectionArea(
-        "POLYGON ((100 300, 305 299, 150 200, 300 150, 150 100, 300 50, 100 50, 100 300))",
-        "POLYGON ((400 350, 150 250, 350 200, 200 150, 350 100, 180 50, 400 50, 400 350))");
-  }
-
   public void testAOverlapBWithHole() {
     checkIntersectionArea(
         "POLYGON ((100 300, 305 299, 150 200, 300 150, 150 100, 300 50, 100 50, 100 300))",
@@ -87,16 +45,8 @@ public void testAOverlapBMultiHole() {
         "MULTIPOLYGON (((55 266, 150 150, 170 290, 55 266)), ((100 0, 70 130, 260 160, 291 45, 100 0), (150 40, 125 98, 220 110, 150 40)))");
   }
 
-  private void checkIntersectionArea(String wktA, String wktB) {
-    Geometry a = read(wktA);
-    Geometry b = read(wktB);
-    
-    OverlayArea ova = new OverlayArea(a);
-    double ovIntArea = ova.intersectionArea(b);
-    
-    double intAreaFull = a.intersection(b).getArea();
-    
-    //System.out.printf("OverlayArea: %f   Full overlay: %f\n", ovIntArea, intAreaFull);
-    assertEquals(intAreaFull, ovIntArea, 0.0001);
+  @Override
+  protected double computeOverlayArea(Geometry a, Geometry b) {
+    return new OverlayArea(a).intersectionArea(b);
   }
 }
diff --git a/modules/lab/src/test/java/org/locationtech/jts/operation/overlayarea/SimpleOverlayAreaTest.java b/modules/lab/src/test/java/org/locationtech/jts/operation/overlayarea/SimpleOverlayAreaTest.java
index c2ab1813d7..bc31d9ad58 100644
--- a/modules/lab/src/test/java/org/locationtech/jts/operation/overlayarea/SimpleOverlayAreaTest.java
+++ b/modules/lab/src/test/java/org/locationtech/jts/operation/overlayarea/SimpleOverlayAreaTest.java
@@ -11,12 +11,11 @@
  */
 package org.locationtech.jts.operation.overlayarea;
 
-import org.locationtech.jts.geom.Polygon;
-
 import junit.textui.TestRunner;
-import test.jts.GeometryTestCase;
+import org.locationtech.jts.geom.Geometry;
+import org.locationtech.jts.geom.Polygon;
 
-public class SimpleOverlayAreaTest extends GeometryTestCase {
+public class SimpleOverlayAreaTest extends BaseOverlayAreaTest {
 
   public static void main(String args[]) {
     TestRunner.run(SimpleOverlayAreaTest.class);
@@ -26,58 +25,8 @@ public SimpleOverlayAreaTest(String name) {
     super(name);
   }
 
-  public void testDisjoint() {
-    checkIntersectionArea(
-        "POLYGON ((10 90, 40 90, 40 60, 10 60, 10 90))",
-        "POLYGON ((90 10, 50 10, 50 50, 90 50, 90 10))");
-  }
-
-  //TODO: fix this bug
-  public void xtestTouching() {
-    checkIntersectionArea(
-        "POLYGON ((10 90, 50 90, 50 50, 10 50, 10 90))",
-        "POLYGON ((90 10, 50 10, 50 50, 90 50, 90 10))");
-  }
-
-  public void testRectangleAContainsB() {
-    checkIntersectionArea(
-        "POLYGON ((100 300, 300 300, 300 100, 100 100, 100 300))",
-        "POLYGON ((150 250, 250 250, 250 150, 150 150, 150 250))");
-  }
-
-  public void testTriangleAContainsB() {
-    checkIntersectionArea(
-        "POLYGON ((60 170, 270 370, 380 60, 60 170))",
-        "POLYGON ((200 250, 245 155, 291 195, 200 250))");
-  }
-
-  public void testRectangleOverlap() {
-    checkIntersectionArea(
-        "POLYGON ((100 200, 200 200, 200 100, 100 100, 100 200))",
-        "POLYGON ((250 250, 250 150, 150 150, 150 250, 250 250))");
-  }
-
-  public void testRectangleTriangleOverlap() {
-    checkIntersectionArea(
-        "POLYGON ((100 200, 200 200, 200 100, 100 100, 100 200))",
-        "POLYGON ((300 200, 150 150, 300 100, 300 200))");
-  }
-
-  public void testSawOverlap() {
-    checkIntersectionArea(
-        "POLYGON ((100 300, 305 299, 150 200, 300 150, 150 100, 300 50, 100 50, 100 300))",
-        "POLYGON ((400 350, 150 250, 350 200, 200 150, 350 100, 180 50, 400 50, 400 350))");
-  }
-
-  private void checkIntersectionArea(String wktA, String wktB) {
-    Polygon a = (Polygon) read(wktA);
-    Polygon b = (Polygon) read(wktB);
-    
-    double ovIntArea = SimpleOverlayArea.intersectionArea(a, b);
-    
-    double intAreaFull = a.intersection(b).getArea();
-    
-    //System.out.printf("OverlayArea: %f   Full overlay: %f\n", ovIntArea, intAreaFull);
-    assertEquals(intAreaFull, ovIntArea, 0.0001);
+  @Override
+  protected double computeOverlayArea(Geometry a, Geometry b) {
+    return SimpleOverlayArea.intersectionArea((Polygon) a, (Polygon) b);
   }
 }

From d918c7e28e7500edeb9613cd2bad138a84f7f25e Mon Sep 17 00:00:00 2001
From: tanelk <tanel.kiis@gmail.com>
Date: Sun, 17 Mar 2024 17:29:14 +0200
Subject: [PATCH 02/12] Symetric tests

Signed-off-by: tanelk <tanel.kiis@gmail.com>
---
 .../operation/overlayarea/OverlayArea.java    |  4 +---
 .../overlayarea/BaseOverlayAreaTest.java      | 21 ++++++++++++-------
 2 files changed, 14 insertions(+), 11 deletions(-)

diff --git a/modules/lab/src/main/java/org/locationtech/jts/operation/overlayarea/OverlayArea.java b/modules/lab/src/main/java/org/locationtech/jts/operation/overlayarea/OverlayArea.java
index 58758ead8f..30766cf253 100644
--- a/modules/lab/src/main/java/org/locationtech/jts/operation/overlayarea/OverlayArea.java
+++ b/modules/lab/src/main/java/org/locationtech/jts/operation/overlayarea/OverlayArea.java
@@ -23,11 +23,9 @@
 import org.locationtech.jts.geom.CoordinateSequence;
 import org.locationtech.jts.geom.Envelope;
 import org.locationtech.jts.geom.Geometry;
-import org.locationtech.jts.geom.GeometryFactory;
 import org.locationtech.jts.geom.GeometryFilter;
 import org.locationtech.jts.geom.LinearRing;
 import org.locationtech.jts.geom.Location;
-import org.locationtech.jts.geom.Point;
 import org.locationtech.jts.geom.Polygon;
 import org.locationtech.jts.index.kdtree.KdNode;
 import org.locationtech.jts.index.kdtree.KdTree;
@@ -179,7 +177,7 @@ private double areaContainedOrDisjoint(LinearRing geom) {
     if (area0 != 0.0) return area0;
     
     // only checking one point, so non-indexed is faster
-    SimplePointInAreaLocator locator = new SimplePointInAreaLocator(geom);
+    SimplePointInAreaLocator locator = new SimplePointInAreaLocator(geom.getFactory().createPolygon(geom));
     double area1 = areaForContainedGeom(geom0, geom.getEnvelopeInternal(), locator);
     // geom0 is either disjoint or contained - either way we are done
     return area1;
diff --git a/modules/lab/src/test/java/org/locationtech/jts/operation/overlayarea/BaseOverlayAreaTest.java b/modules/lab/src/test/java/org/locationtech/jts/operation/overlayarea/BaseOverlayAreaTest.java
index 01feff1ea0..68c94c7a14 100644
--- a/modules/lab/src/test/java/org/locationtech/jts/operation/overlayarea/BaseOverlayAreaTest.java
+++ b/modules/lab/src/test/java/org/locationtech/jts/operation/overlayarea/BaseOverlayAreaTest.java
@@ -27,48 +27,48 @@ public BaseOverlayAreaTest(String name) {
   }
 
   public void testDisjoint() {
-    checkIntersectionArea(
+    checkIntersectionAreaSymmetric(
         "POLYGON ((10 90, 40 90, 40 60, 10 60, 10 90))",
         "POLYGON ((90 10, 50 10, 50 50, 90 50, 90 10))");
   }
 
   public void testTouching() {
-    checkIntersectionArea(
+    checkIntersectionAreaSymmetric(
         "POLYGON ((10 90, 50 90, 50 50, 10 50, 10 90))",
         "POLYGON ((90 10, 50 10, 50 50, 90 50, 90 10))");
   }
 
   public void testRectangleAContainsB() {
-    checkIntersectionArea(
+    checkIntersectionAreaSymmetric(
         "POLYGON ((100 300, 300 300, 300 100, 100 100, 100 300))",
         "POLYGON ((150 250, 250 250, 250 150, 150 150, 150 250))");
   }
 
   public void testTriangleAContainsB() {
-    checkIntersectionArea(
+    checkIntersectionAreaSymmetric(
         "POLYGON ((60 170, 270 370, 380 60, 60 170))",
         "POLYGON ((200 250, 245 155, 291 195, 200 250))");
   }
 
   public void testRectangleOverlap() {
-    checkIntersectionArea(
+    checkIntersectionAreaSymmetric(
         "POLYGON ((100 200, 200 200, 200 100, 100 100, 100 200))",
         "POLYGON ((250 250, 250 150, 150 150, 150 250, 250 250))");
   }
 
   public void testRectangleTriangleOverlap() {
-    checkIntersectionArea(
+    checkIntersectionAreaSymmetric(
         "POLYGON ((100 200, 200 200, 200 100, 100 100, 100 200))",
         "POLYGON ((300 200, 150 150, 300 100, 300 200))");
   }
 
   public void testSawOverlap() {
-    checkIntersectionArea(
+    checkIntersectionAreaSymmetric(
         "POLYGON ((100 300, 305 299, 150 200, 300 150, 150 100, 300 50, 100 50, 100 300))",
         "POLYGON ((400 350, 150 250, 350 200, 200 150, 350 100, 180 50, 400 50, 400 350))");
   }
 
-  protected void checkIntersectionArea(String wktA, String wktB) {
+  protected final void checkIntersectionArea(String wktA, String wktB) {
     Geometry a = read(wktA);
     Geometry b = read(wktB);
     
@@ -80,5 +80,10 @@ protected void checkIntersectionArea(String wktA, String wktB) {
     assertEquals(intAreaFull, ovIntArea, 0.0001);
   }
 
+  protected final void checkIntersectionAreaSymmetric(String wktA, String wktB) {
+    checkIntersectionArea(wktA, wktB);
+    checkIntersectionArea(wktB, wktA);
+  }
+
   abstract protected double computeOverlayArea(Geometry a, Geometry b);
 }

From 79bbbef78329d0b20d9db94ec7c3e877dc2b1dbc Mon Sep 17 00:00:00 2001
From: tanelk <tanel.kiis@gmail.com>
Date: Sun, 17 Mar 2024 17:54:27 +0200
Subject: [PATCH 03/12] Add more testcases

Signed-off-by: tanelk <tanel.kiis@gmail.com>
---
 .../overlayarea/BaseOverlayAreaTest.java      | 30 +++++++++++++++++++
 1 file changed, 30 insertions(+)

diff --git a/modules/lab/src/test/java/org/locationtech/jts/operation/overlayarea/BaseOverlayAreaTest.java b/modules/lab/src/test/java/org/locationtech/jts/operation/overlayarea/BaseOverlayAreaTest.java
index 68c94c7a14..98c7ecf976 100644
--- a/modules/lab/src/test/java/org/locationtech/jts/operation/overlayarea/BaseOverlayAreaTest.java
+++ b/modules/lab/src/test/java/org/locationtech/jts/operation/overlayarea/BaseOverlayAreaTest.java
@@ -38,6 +38,12 @@ public void testTouching() {
         "POLYGON ((90 10, 50 10, 50 50, 90 50, 90 10))");
   }
 
+  public void testTouchingNonPerpendicular() {
+    checkIntersectionAreaSymmetric(
+            "POLYGON ((10 90, 49 90, 50 50, 10 53, 10 90))",
+            "POLYGON ((90 10, 47 10, 50 50, 90 54, 90 10))");
+  }
+
   public void testRectangleAContainsB() {
     checkIntersectionAreaSymmetric(
         "POLYGON ((100 300, 300 300, 300 100, 100 100, 100 300))",
@@ -68,6 +74,30 @@ public void testSawOverlap() {
         "POLYGON ((400 350, 150 250, 350 200, 200 150, 350 100, 180 50, 400 50, 400 350))");
   }
 
+  public void testVertexIntersectionContained() {
+    checkIntersectionAreaSymmetric(
+        "POLYGON ((10 10, 40 50, 90 10, 10 10))",
+            "POLYGON ((30 20, 40 50, 60 20, 30 20))");
+  }
+
+  public void testVertexIntersectionOverlapping() {
+    checkIntersectionAreaSymmetric(
+            "POLYGON ((10 10, 40 50, 90 10, 10 10))",
+            "POLYGON ((30 -20, 40 50, 60 -20, 30 -20))");
+  }
+
+  public void testVertexIntersectionOverlapping2() {
+    checkIntersectionAreaSymmetric(
+            "POLYGON ((10 10, 40 50, 90 10, 10 10))",
+            "POLYGON ((30 20, 40 50, 90 20, 30 20))");
+  }
+
+  public void testVertexIntersectionTwoAreas() {
+    checkIntersectionAreaSymmetric(
+            "POLYGON ((10 10, 30 50, 60 10, 10 10))",
+            "POLYGON ((20 60, 40 60, 40 20, 30 50, 20 20, 20 60))");
+  }
+
   protected final void checkIntersectionArea(String wktA, String wktB) {
     Geometry a = read(wktA);
     Geometry b = read(wktB);

From 346db8902015d3b46e64d54f907830ec277b5f25 Mon Sep 17 00:00:00 2001
From: tanelk <tanel.kiis@gmail.com>
Date: Sun, 17 Mar 2024 18:25:43 +0200
Subject: [PATCH 04/12] Add more testcases

Signed-off-by: tanelk <tanel.kiis@gmail.com>
---
 .../operation/overlayarea/OverlayArea.java    |  7 +++--
 .../overlayarea/BaseOverlayAreaTest.java      | 26 ++++++++++++++++---
 2 files changed, 27 insertions(+), 6 deletions(-)

diff --git a/modules/lab/src/main/java/org/locationtech/jts/operation/overlayarea/OverlayArea.java b/modules/lab/src/main/java/org/locationtech/jts/operation/overlayarea/OverlayArea.java
index 30766cf253..e3fda0ac00 100644
--- a/modules/lab/src/main/java/org/locationtech/jts/operation/overlayarea/OverlayArea.java
+++ b/modules/lab/src/main/java/org/locationtech/jts/operation/overlayarea/OverlayArea.java
@@ -275,12 +275,15 @@ public void processIntersections(SegmentString a, int aIndex, SegmentString b, i
         }
 
       } else {
-        // vertex-vertex intersection
+        // vertex-vertex inteqrsection
         // This intersection is visited 4 times - include only once
-        if (!intPt.equals2D(a1) || !intPt.equals2D(b1)) {
+        if (!a1.equals2D(b1)) {
           return;
         }
 
+        // If A0->A1 is collinear with B0->B1, then the intersection point might not be equal to A1 and B1
+        intPt = a1;
+
         Coordinate a2 = a.nextInRing(aIndex + 1);
         Coordinate b2 = b.nextInRing(bIndex + 1);
         if (isCCWA) {
diff --git a/modules/lab/src/test/java/org/locationtech/jts/operation/overlayarea/BaseOverlayAreaTest.java b/modules/lab/src/test/java/org/locationtech/jts/operation/overlayarea/BaseOverlayAreaTest.java
index 98c7ecf976..ff60b7bc5c 100644
--- a/modules/lab/src/test/java/org/locationtech/jts/operation/overlayarea/BaseOverlayAreaTest.java
+++ b/modules/lab/src/test/java/org/locationtech/jts/operation/overlayarea/BaseOverlayAreaTest.java
@@ -80,6 +80,24 @@ public void testVertexIntersectionContained() {
             "POLYGON ((30 20, 40 50, 60 20, 30 20))");
   }
 
+  public void testVertexIntersectionContained2() {
+    checkIntersectionAreaSymmetric(
+            "POLYGON ((10 10, 40 50, 90 10, 10 10))",
+            "POLYGON ((30 10, 40 50, 60 10, 30 10))");
+  }
+
+  public void testVertexIntersectionContained3() {
+    checkIntersectionAreaSymmetric(
+            "POLYGON ((10 10, 40 50, 90 10, 10 10))",
+            "POLYGON ((10 10, 40 50, 60 20, 10 10))");
+  }
+
+  public void testVertexIntersectionContained4() {
+    checkIntersectionAreaSymmetric(
+            "POLYGON ((10 10, 10 50, 50 50, 50 10, 10 10))",
+            "POLYGON ((10 10, 30 30, 50 10, 10 10))");
+  }
+
   public void testVertexIntersectionOverlapping() {
     checkIntersectionAreaSymmetric(
             "POLYGON ((10 10, 40 50, 90 10, 10 10))",
@@ -102,12 +120,12 @@ protected final void checkIntersectionArea(String wktA, String wktB) {
     Geometry a = read(wktA);
     Geometry b = read(wktB);
     
-    double ovIntArea = computeOverlayArea(a, b);
-    
     double intAreaFull = a.intersection(b).getArea();
     
-    //System.out.printf("OverlayArea: %f   Full overlay: %f\n", ovIntArea, intAreaFull);
-    assertEquals(intAreaFull, ovIntArea, 0.0001);
+    assertEquals(intAreaFull, computeOverlayArea(a, b), 0.0001);
+    assertEquals(intAreaFull, computeOverlayArea(a.reverse(), b), 0.0001);
+    assertEquals(intAreaFull, computeOverlayArea(a, b.reverse()), 0.0001);
+    assertEquals(intAreaFull, computeOverlayArea(a.reverse(), b.reverse()), 0.0001);
   }
 
   protected final void checkIntersectionAreaSymmetric(String wktA, String wktB) {

From 2543ce3137f793f2feffde11d5d0949354399e38 Mon Sep 17 00:00:00 2001
From: tanelk <tanel.kiis@gmail.com>
Date: Sun, 17 Mar 2024 18:57:02 +0200
Subject: [PATCH 05/12] Fix typo

Signed-off-by: tanelk <tanel.kiis@gmail.com>
---
 .../locationtech/jts/operation/overlayarea/OverlayArea.java | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/modules/lab/src/main/java/org/locationtech/jts/operation/overlayarea/OverlayArea.java b/modules/lab/src/main/java/org/locationtech/jts/operation/overlayarea/OverlayArea.java
index e3fda0ac00..91b5e08d60 100644
--- a/modules/lab/src/main/java/org/locationtech/jts/operation/overlayarea/OverlayArea.java
+++ b/modules/lab/src/main/java/org/locationtech/jts/operation/overlayarea/OverlayArea.java
@@ -275,7 +275,7 @@ public void processIntersections(SegmentString a, int aIndex, SegmentString b, i
         }
 
       } else {
-        // vertex-vertex inteqrsection
+        // vertex-vertex intersection
         // This intersection is visited 4 times - include only once
         if (!a1.equals2D(b1)) {
           return;
@@ -296,8 +296,8 @@ public void processIntersections(SegmentString a, int aIndex, SegmentString b, i
         double aAngle = Angle.interiorAngle(a0, intPt, a2);
         double bAngle = Angle.interiorAngle(b0, intPt, b2);
 
-        // The LTE ja LT are chosen such that when A0->A1 is collinear with B0->B1, then A is chosen
-        // and when A1->A2 is collinear with B1->B2, then A is chosen.
+        // The LTE ja LT are chosen such that when A0->A1 is collinear with B0->B1,
+        // or when A1->A2 is collinear with B1->B2, then A is chosen.
         // This avoids double counting in the case of collinear segments.
         if (Angle.interiorAngle(a0, intPt, b2) <= bAngle) {
           area += EdgeVector.area2Term(intPt, a1, a0, false);

From cf63d7e7376f70b91b365ffc68c147259f0c4c2b Mon Sep 17 00:00:00 2001
From: tanelk <tanel.kiis@gmail.com>
Date: Sun, 17 Mar 2024 20:41:21 +0200
Subject: [PATCH 06/12] Fix vertex-edge case

Signed-off-by: tanelk <tanel.kiis@gmail.com>
---
 .../jts/operation/overlayarea/OverlayArea.java       | 10 ++++++++++
 .../operation/overlayarea/BaseOverlayAreaTest.java   | 12 ++++++++++++
 2 files changed, 22 insertions(+)

diff --git a/modules/lab/src/main/java/org/locationtech/jts/operation/overlayarea/OverlayArea.java b/modules/lab/src/main/java/org/locationtech/jts/operation/overlayarea/OverlayArea.java
index 91b5e08d60..e5ce206ee2 100644
--- a/modules/lab/src/main/java/org/locationtech/jts/operation/overlayarea/OverlayArea.java
+++ b/modules/lab/src/main/java/org/locationtech/jts/operation/overlayarea/OverlayArea.java
@@ -267,9 +267,19 @@ public void processIntersections(SegmentString a, int aIndex, SegmentString b, i
          */
 
         if (Orientation.CLOCKWISE == Orientation.index(a0, a1, b0)) {
+          if (intPt.equals2D(a1)) {
+            // Intersection at vertex and A0 -> A1 is outside the intersection area.
+            // Area will be computed by the segment A1 -> A2
+            return;
+          }
           area += EdgeVector.area2Term(intPt, a0, a1, true);
           area += EdgeVector.area2Term(intPt, b1, b0, false);
         } else if (Orientation.CLOCKWISE == Orientation.index(a0, a1, b1)) {
+          if (intPt.equals2D(a0)) {
+            // Intersection at vertex and A0 -> A1 is outside the intersection area.
+            // Area will be computed by the segment A(-1) -> A0
+            return;
+          }
           area += EdgeVector.area2Term(intPt, a1, a0, false);
           area += EdgeVector.area2Term(intPt, b0, b1, true);
         }
diff --git a/modules/lab/src/test/java/org/locationtech/jts/operation/overlayarea/BaseOverlayAreaTest.java b/modules/lab/src/test/java/org/locationtech/jts/operation/overlayarea/BaseOverlayAreaTest.java
index ff60b7bc5c..fb21729a99 100644
--- a/modules/lab/src/test/java/org/locationtech/jts/operation/overlayarea/BaseOverlayAreaTest.java
+++ b/modules/lab/src/test/java/org/locationtech/jts/operation/overlayarea/BaseOverlayAreaTest.java
@@ -116,6 +116,18 @@ public void testVertexIntersectionTwoAreas() {
             "POLYGON ((20 60, 40 60, 40 20, 30 50, 20 20, 20 60))");
   }
 
+  public void testVertexIntersectionOnEdge() {
+    checkIntersectionAreaSymmetric(
+            "POLYGON ((10 10, 10 40, 40 40, 40 10, 10 10))",
+            "POLYGON ((20 30, 20 40, 70 50, 20 30))");
+  }
+
+  public void testVertexIntersectionOnEdge2() {
+    checkIntersectionAreaSymmetric(
+            "POLYGON ((10 30, 10 60, 40 60, 40 30, 10 30))",
+            "POLYGON ((40 10, 20 30, 40 50, 50 30, 40 10))");
+  }
+
   protected final void checkIntersectionArea(String wktA, String wktB) {
     Geometry a = read(wktA);
     Geometry b = read(wktB);

From 1d70ba65ba06ea5a31ce5bde808342bc81455aac Mon Sep 17 00:00:00 2001
From: tanelk <tanel.kiis@gmail.com>
Date: Mon, 18 Mar 2024 07:28:45 +0200
Subject: [PATCH 07/12] Better comments

Signed-off-by: tanelk <tanel.kiis@gmail.com>
---
 .../operation/overlayarea/OverlayArea.java    | 22 +++++++++++++++----
 1 file changed, 18 insertions(+), 4 deletions(-)

diff --git a/modules/lab/src/main/java/org/locationtech/jts/operation/overlayarea/OverlayArea.java b/modules/lab/src/main/java/org/locationtech/jts/operation/overlayarea/OverlayArea.java
index e5ce206ee2..391ae4aa06 100644
--- a/modules/lab/src/main/java/org/locationtech/jts/operation/overlayarea/OverlayArea.java
+++ b/modules/lab/src/main/java/org/locationtech/jts/operation/overlayarea/OverlayArea.java
@@ -291,9 +291,14 @@ public void processIntersections(SegmentString a, int aIndex, SegmentString b, i
           return;
         }
 
-        // If A0->A1 is collinear with B0->B1, then the intersection point might not be equal to A1 and B1
+        // If A0->A1 is collinear with B0->B1,
+        // then the intersection point from LineIntersector might not be equal to A1 and B1
         intPt = a1;
 
+        /* Get the next vertices in the CW direction.
+        Now we have four segments: A0->A1, A1->A2, B0->B1, B1->B2
+        and the intersection point is A1 == B1.
+         */
         Coordinate a2 = a.nextInRing(aIndex + 1);
         Coordinate b2 = b.nextInRing(bIndex + 1);
         if (isCCWA) {
@@ -303,12 +308,21 @@ public void processIntersections(SegmentString a, int aIndex, SegmentString b, i
           b2 = b.prevInRing(bIndex);
         }
 
+        /* The angles A0->A1->A2 and B0->B1->B2 determine
+         the maximum intersection area interior angle.
+         Edges from the other polygon that lie within this angle
+         are on the boundary of the intersection area.
+
+         Depending on the relative orientation of the polygons,
+         we could pick 0, 2 or 4 segments to contribute to the area.
+
+        The LTE ja LT are chosen such that when A0->A1 is collinear with B0->B1,
+        or when A1->A2 is collinear with B1->B2, then only the segment from polygon A
+        is chosen to avoid double counting.
+         */
         double aAngle = Angle.interiorAngle(a0, intPt, a2);
         double bAngle = Angle.interiorAngle(b0, intPt, b2);
 
-        // The LTE ja LT are chosen such that when A0->A1 is collinear with B0->B1,
-        // or when A1->A2 is collinear with B1->B2, then A is chosen.
-        // This avoids double counting in the case of collinear segments.
         if (Angle.interiorAngle(a0, intPt, b2) <= bAngle) {
           area += EdgeVector.area2Term(intPt, a1, a0, false);
         }

From ca380fe6a00760b564b0481491fe0eac32e3cba4 Mon Sep 17 00:00:00 2001
From: tanelk <tanel.kiis@gmail.com>
Date: Mon, 18 Mar 2024 08:33:36 +0200
Subject: [PATCH 08/12] Avoid recalc

Signed-off-by: tanelk <tanel.kiis@gmail.com>
---
 .../jts/operation/overlayarea/OverlayArea.java    | 15 +++++++++------
 1 file changed, 9 insertions(+), 6 deletions(-)

diff --git a/modules/lab/src/main/java/org/locationtech/jts/operation/overlayarea/OverlayArea.java b/modules/lab/src/main/java/org/locationtech/jts/operation/overlayarea/OverlayArea.java
index 391ae4aa06..12b71084d9 100644
--- a/modules/lab/src/main/java/org/locationtech/jts/operation/overlayarea/OverlayArea.java
+++ b/modules/lab/src/main/java/org/locationtech/jts/operation/overlayarea/OverlayArea.java
@@ -320,19 +320,22 @@ public void processIntersections(SegmentString a, int aIndex, SegmentString b, i
         or when A1->A2 is collinear with B1->B2, then only the segment from polygon A
         is chosen to avoid double counting.
          */
-        double aAngle = Angle.interiorAngle(a0, intPt, a2);
-        double bAngle = Angle.interiorAngle(b0, intPt, b2);
+        double aaAngle = Angle.interiorAngle(a0, intPt, a2);
+        double bbAngle = Angle.interiorAngle(b0, intPt, b2);
 
-        if (Angle.interiorAngle(a0, intPt, b2) <= bAngle) {
+        double abAngle = Angle.interiorAngle(a0, intPt, b2);
+        double baAngle = Angle.interiorAngle(b0, intPt, a2);
+
+        if (abAngle <= bbAngle) {
           area += EdgeVector.area2Term(intPt, a1, a0, false);
         }
-        if (Angle.interiorAngle(b0, intPt, a2) <= bAngle) {
+        if (baAngle <= bbAngle) {
           area += EdgeVector.area2Term(intPt, a1, a2, true);
         }
-        if (Angle.interiorAngle(b0, intPt, a2) < aAngle) {
+        if (baAngle < aaAngle) {
           area += EdgeVector.area2Term(intPt, b1, b0, false);
         }
-        if (Angle.interiorAngle(a0, intPt, b2) < aAngle) {
+        if (abAngle < aaAngle) {
           area += EdgeVector.area2Term(intPt, b1, b2, true);
         }
       }

From 1df695e0e19126666029f307b7aa1f1136c53441 Mon Sep 17 00:00:00 2001
From: tanelk <tanel.kiis@gmail.com>
Date: Mon, 18 Mar 2024 08:51:39 +0200
Subject: [PATCH 09/12] Better naming

Signed-off-by: tanelk <tanel.kiis@gmail.com>
---
 .../jts/operation/overlayarea/OverlayArea.java   | 16 ++++++++--------
 1 file changed, 8 insertions(+), 8 deletions(-)

diff --git a/modules/lab/src/main/java/org/locationtech/jts/operation/overlayarea/OverlayArea.java b/modules/lab/src/main/java/org/locationtech/jts/operation/overlayarea/OverlayArea.java
index 12b71084d9..73ae602f07 100644
--- a/modules/lab/src/main/java/org/locationtech/jts/operation/overlayarea/OverlayArea.java
+++ b/modules/lab/src/main/java/org/locationtech/jts/operation/overlayarea/OverlayArea.java
@@ -320,22 +320,22 @@ public void processIntersections(SegmentString a, int aIndex, SegmentString b, i
         or when A1->A2 is collinear with B1->B2, then only the segment from polygon A
         is chosen to avoid double counting.
          */
-        double aaAngle = Angle.interiorAngle(a0, intPt, a2);
-        double bbAngle = Angle.interiorAngle(b0, intPt, b2);
+        double angleA0A2 = Angle.interiorAngle(a0, intPt, a2);
+        double angleB0B2 = Angle.interiorAngle(b0, intPt, b2);
 
-        double abAngle = Angle.interiorAngle(a0, intPt, b2);
-        double baAngle = Angle.interiorAngle(b0, intPt, a2);
+        double angleA0B2 = Angle.interiorAngle(a0, intPt, b2);
+        double angleB0A2 = Angle.interiorAngle(b0, intPt, a2);
 
-        if (abAngle <= bbAngle) {
+        if (angleA0B2 <= angleB0B2) {
           area += EdgeVector.area2Term(intPt, a1, a0, false);
         }
-        if (baAngle <= bbAngle) {
+        if (angleB0A2 <= angleB0B2) {
           area += EdgeVector.area2Term(intPt, a1, a2, true);
         }
-        if (baAngle < aaAngle) {
+        if (angleB0A2 < angleA0A2) {
           area += EdgeVector.area2Term(intPt, b1, b0, false);
         }
-        if (abAngle < aaAngle) {
+        if (angleA0B2 < angleA0A2) {
           area += EdgeVector.area2Term(intPt, b1, b2, true);
         }
       }

From 59cf2461459b064a2c45040e42af4371739e14ec Mon Sep 17 00:00:00 2001
From: tanelk <tanel.kiis@gmail.com>
Date: Mon, 18 Mar 2024 10:02:16 +0200
Subject: [PATCH 10/12] Tests for collinear edges

Signed-off-by: tanelk <tanel.kiis@gmail.com>
---
 .../operation/overlayarea/BaseOverlayAreaTest.java | 14 ++++++++++++++
 1 file changed, 14 insertions(+)

diff --git a/modules/lab/src/test/java/org/locationtech/jts/operation/overlayarea/BaseOverlayAreaTest.java b/modules/lab/src/test/java/org/locationtech/jts/operation/overlayarea/BaseOverlayAreaTest.java
index fb21729a99..17550f30d5 100644
--- a/modules/lab/src/test/java/org/locationtech/jts/operation/overlayarea/BaseOverlayAreaTest.java
+++ b/modules/lab/src/test/java/org/locationtech/jts/operation/overlayarea/BaseOverlayAreaTest.java
@@ -128,6 +128,20 @@ public void testVertexIntersectionOnEdge2() {
             "POLYGON ((40 10, 20 30, 40 50, 50 30, 40 10))");
   }
 
+  public void testCollinearOverlappingEdgesPartial() {
+    checkIntersectionAreaSymmetric(
+        "POLYGON ((10 30, 30 30, 30 10, 10 10, 10 30))",
+        "POLYGON ((20 30, 40 30, 40 10, 20 10, 20 30))"
+    );
+  }
+
+  public void testCollinearOverlappingEdgesFull() {
+    checkIntersectionAreaSymmetric(
+        "POLYGON ((10 30, 50 30, 50 10, 10 10, 10 30))",
+        "POLYGON ((20 30, 40 30, 40 10, 20 10, 20 30))"
+    );
+  }
+
   protected final void checkIntersectionArea(String wktA, String wktB) {
     Geometry a = read(wktA);
     Geometry b = read(wktB);

From 6d2bd0ac3a232df7eb7542ebcb2884a7dcd305e1 Mon Sep 17 00:00:00 2001
From: tanelk <tanel.kiis@gmail.com>
Date: Fri, 22 Mar 2024 21:56:20 +0200
Subject: [PATCH 11/12] Fix SimpleOverlayArea

Signed-off-by: tanelk <tanel.kiis@gmail.com>
---
 .../overlayarea/IntersectionVisitor.java      | 142 ++++++++++++++++++
 .../operation/overlayarea/OverlayArea.java    | 119 +--------------
 .../overlayarea/SimpleOverlayArea.java        |  69 ++-------
 3 files changed, 157 insertions(+), 173 deletions(-)
 create mode 100644 modules/lab/src/main/java/org/locationtech/jts/operation/overlayarea/IntersectionVisitor.java

diff --git a/modules/lab/src/main/java/org/locationtech/jts/operation/overlayarea/IntersectionVisitor.java b/modules/lab/src/main/java/org/locationtech/jts/operation/overlayarea/IntersectionVisitor.java
new file mode 100644
index 0000000000..b81b5b1e4e
--- /dev/null
+++ b/modules/lab/src/main/java/org/locationtech/jts/operation/overlayarea/IntersectionVisitor.java
@@ -0,0 +1,142 @@
+/*
+ * Copyright (c) 2022 Martin Davis, and others.
+ *
+ * All rights reserved. This program and the accompanying materials
+ * are made available under the terms of the Eclipse Public License 2.0
+ * and Eclipse Distribution License v. 1.0 which accompanies this distribution.
+ * The Eclipse Public License is available at http://www.eclipse.org/legal/epl-v20.html
+ * and the Eclipse Distribution License is available at
+ *
+ * http://www.eclipse.org/org/documents/edl-v10.php.
+ */
+package org.locationtech.jts.operation.overlayarea;
+
+import org.locationtech.jts.algorithm.Angle;
+import org.locationtech.jts.algorithm.LineIntersector;
+import org.locationtech.jts.algorithm.Orientation;
+import org.locationtech.jts.algorithm.RobustLineIntersector;
+import org.locationtech.jts.geom.Coordinate;
+import org.locationtech.jts.noding.SegmentIntersector;
+import org.locationtech.jts.noding.SegmentString;
+
+/**
+ * Computes the partial overlay area of two polygons by summing the contributions of
+ * the edge vectors created by the intersections between the edges of the two polygons.
+ */
+class IntersectionVisitor implements SegmentIntersector {
+
+    private static final LineIntersector li = new RobustLineIntersector();
+
+    private double area = 0.0;
+
+    double getArea() {
+        return area;
+    }
+
+    @Override
+    public void processIntersections(SegmentString a, int aIndex, SegmentString b, int bIndex) {
+        boolean isCCWA = (boolean) a.getData();
+        boolean isCCWB = (boolean) b.getData();
+
+        Coordinate a0 = a.getCoordinate(aIndex);
+        Coordinate a1 = a.getCoordinate(aIndex + 1);
+        Coordinate b0 = b.getCoordinate(bIndex);
+        Coordinate b1 = b.getCoordinate(bIndex + 1);
+
+        if (isCCWA) {
+            Coordinate tmp = a0;
+            a0 = a1;
+            a1 = tmp;
+        }
+        if (isCCWB) {
+            Coordinate tmp = b0;
+            b0 = b1;
+            b1 = tmp;
+        }
+
+        li.computeIntersection(a0, a1, b0, b1);
+        if (!li.hasIntersection()) return;
+
+        if (li.isProper() || li.isInteriorIntersection()) {
+            // Edge-edge intersection OR vertex-edge intersection
+
+            /**
+             * An intersection creates two edge vectors which contribute to the area.
+             *
+             * With both rings oriented CW (effectively)
+             * There are two situations for segment intersection:
+             *
+             * 1) A entering B, B exiting A => rays are IP->A1:R, IP->B0:L
+             * 2) A exiting B, B entering A => rays are IP->A0:L, IP->B1:R
+             * (where IP is the intersection point,
+             * and  :L/R indicates result polygon interior is to the Left or Right).
+             *
+             * For accuracy the full edge is used to provide the direction vector.
+             */
+
+            Coordinate intPt = li.getIntersection(0);
+
+            if (Orientation.CLOCKWISE == Orientation.index(a0, a1, b0)) {
+                if (intPt.equals2D(a1)) {
+                    // Intersection at vertex and A0 -> A1 is outside the intersection area.
+                    // Area will be computed by the segment A1 -> A2
+                    return;
+                }
+                area += EdgeVector.area2Term(intPt, a0, a1, true);
+                area += EdgeVector.area2Term(intPt, b1, b0, false);
+            } else if (Orientation.CLOCKWISE == Orientation.index(a0, a1, b1)) {
+                if (intPt.equals2D(a0)) {
+                    // Intersection at vertex and A0 -> A1 is outside the intersection area.
+                    // Area will be computed by the segment A(-1) -> A0
+                    return;
+                }
+                area += EdgeVector.area2Term(intPt, a1, a0, false);
+                area += EdgeVector.area2Term(intPt, b0, b1, true);
+            }
+
+        } else {
+            // vertex-vertex intersection
+            // This intersection is visited 4 times - include only once
+            if (!a1.equals2D(b1)) {
+                return;
+            }
+
+            // If A0->A1 is collinear with B0->B1, then the intersection point from LineIntersector might not be equal to A1 and B1
+            Coordinate intPt = a1;
+
+            Coordinate a2 = a.nextInRing(aIndex + 1);
+            Coordinate b2 = b.nextInRing(bIndex + 1);
+            if (isCCWA) {
+                a2 = a.prevInRing(aIndex);
+            }
+            if (isCCWB) {
+                b2 = b.prevInRing(bIndex);
+            }
+
+            double aAngle = Angle.interiorAngle(a0, intPt, a2);
+            double bAngle = Angle.interiorAngle(b0, intPt, b2);
+
+            // The LTE ja LT are chosen such that when A0->A1 is collinear with B0->B1,
+            // or when A1->A2 is collinear with B1->B2, then A is chosen.
+            // This avoids double counting in the case of collinear segments.
+            if (Angle.interiorAngle(a0, intPt, b2) <= bAngle) {
+                area += EdgeVector.area2Term(intPt, a1, a0, false);
+            }
+            if (Angle.interiorAngle(b0, intPt, a2) <= bAngle) {
+                area += EdgeVector.area2Term(intPt, a1, a2, true);
+            }
+            if (Angle.interiorAngle(b0, intPt, a2) < aAngle) {
+                area += EdgeVector.area2Term(intPt, b1, b0, false);
+            }
+            if (Angle.interiorAngle(a0, intPt, b2) < aAngle) {
+                area += EdgeVector.area2Term(intPt, b1, b2, true);
+            }
+        }
+    }
+
+    @Override
+    public boolean isDone() {
+        // Process all intersections
+        return false;
+    }
+}
diff --git a/modules/lab/src/main/java/org/locationtech/jts/operation/overlayarea/OverlayArea.java b/modules/lab/src/main/java/org/locationtech/jts/operation/overlayarea/OverlayArea.java
index e5ce206ee2..defb781d97 100644
--- a/modules/lab/src/main/java/org/locationtech/jts/operation/overlayarea/OverlayArea.java
+++ b/modules/lab/src/main/java/org/locationtech/jts/operation/overlayarea/OverlayArea.java
@@ -11,11 +11,8 @@
  */
 package org.locationtech.jts.operation.overlayarea;
 
-import org.locationtech.jts.algorithm.Angle;
 import org.locationtech.jts.algorithm.Area;
-import org.locationtech.jts.algorithm.LineIntersector;
 import org.locationtech.jts.algorithm.Orientation;
-import org.locationtech.jts.algorithm.RobustLineIntersector;
 import org.locationtech.jts.algorithm.locate.IndexedPointInAreaLocator;
 import org.locationtech.jts.algorithm.locate.PointOnGeometryLocator;
 import org.locationtech.jts.algorithm.locate.SimplePointInAreaLocator;
@@ -32,7 +29,6 @@
 import org.locationtech.jts.math.MathUtil;
 import org.locationtech.jts.noding.BasicSegmentString;
 import org.locationtech.jts.noding.MCIndexSegmentSetMutualIntersector;
-import org.locationtech.jts.noding.SegmentIntersector;
 import org.locationtech.jts.noding.SegmentSetMutualIntersector;
 import org.locationtech.jts.noding.SegmentString;
 
@@ -43,7 +39,7 @@
  * Computes the area of the overlay of two polygons without forming
  * the actual topology of the overlay.
  * Since the topology is not needed, the computation is
- * is insensitive to the fine details of the overlay topology,
+ * insensitive to the fine details of the overlay topology,
  * and hence is fully robust.
  * It also allows for a simpler implementation with more aggressive
  * performance optimization.
@@ -78,8 +74,6 @@ private static boolean interacts(Geometry geom0, Geometry geom1) {
     return geom0.getEnvelopeInternal().intersects(geom1.getEnvelopeInternal());
   }
 
-  private static LineIntersector li = new RobustLineIntersector();
-  
   private Geometry geom0;
   private Envelope geomEnv0;
   private IndexedPointInAreaLocator locator0;
@@ -220,117 +214,6 @@ private double areaForIntersections(LinearRing geom) {
     return intVisitor.getArea();
   }
 
-  class IntersectionVisitor implements SegmentIntersector {
-    private double area = 0.0;
-
-    double getArea() {
-      return area;
-    }
-
-    @Override
-    public void processIntersections(SegmentString a, int aIndex, SegmentString b, int bIndex) {
-      boolean isCCWA = (boolean) a.getData();
-      boolean isCCWB = (boolean) b.getData();
-
-      Coordinate a0 = a.getCoordinate(aIndex);
-      Coordinate a1 = a.getCoordinate(aIndex + 1);
-      Coordinate b0 = b.getCoordinate(bIndex);
-      Coordinate b1 = b.getCoordinate(bIndex + 1);
-
-      if (isCCWA) {
-        Coordinate tmp = a0; a0 = a1; a1 = tmp;
-      }
-      if (isCCWB) {
-          Coordinate tmp = b0; b0 = b1; b1 = tmp;
-      }
-
-      li.computeIntersection(a0, a1, b0, b1);
-      if (! li.hasIntersection()) return;
-
-      Coordinate intPt = li.getIntersection(0);
-
-      if (li.isProper() || li.isInteriorIntersection()) {
-        // Edge-edge intersection OR vertex-edge intersection
-
-        /**
-         * An intersection creates two edge vectors which contribute to the area.
-         *
-         * With both rings oriented CW (effectively)
-         * There are two situations for segment intersection:
-         *
-         * 1) A entering B, B exiting A => rays are IP->A1:R, IP->B0:L
-         * 2) A exiting B, B entering A => rays are IP->A0:L, IP->B1:R
-         * (where IP is the intersection point,
-         * and  :L/R indicates result polygon interior is to the Left or Right).
-         *
-         * For accuracy the full edge is used to provide the direction vector.
-         */
-
-        if (Orientation.CLOCKWISE == Orientation.index(a0, a1, b0)) {
-          if (intPt.equals2D(a1)) {
-            // Intersection at vertex and A0 -> A1 is outside the intersection area.
-            // Area will be computed by the segment A1 -> A2
-            return;
-          }
-          area += EdgeVector.area2Term(intPt, a0, a1, true);
-          area += EdgeVector.area2Term(intPt, b1, b0, false);
-        } else if (Orientation.CLOCKWISE == Orientation.index(a0, a1, b1)) {
-          if (intPt.equals2D(a0)) {
-            // Intersection at vertex and A0 -> A1 is outside the intersection area.
-            // Area will be computed by the segment A(-1) -> A0
-            return;
-          }
-          area += EdgeVector.area2Term(intPt, a1, a0, false);
-          area += EdgeVector.area2Term(intPt, b0, b1, true);
-        }
-
-      } else {
-        // vertex-vertex intersection
-        // This intersection is visited 4 times - include only once
-        if (!a1.equals2D(b1)) {
-          return;
-        }
-
-        // If A0->A1 is collinear with B0->B1, then the intersection point might not be equal to A1 and B1
-        intPt = a1;
-
-        Coordinate a2 = a.nextInRing(aIndex + 1);
-        Coordinate b2 = b.nextInRing(bIndex + 1);
-        if (isCCWA) {
-          a2 = a.prevInRing(aIndex);
-        }
-        if (isCCWB) {
-          b2 = b.prevInRing(bIndex);
-        }
-
-        double aAngle = Angle.interiorAngle(a0, intPt, a2);
-        double bAngle = Angle.interiorAngle(b0, intPt, b2);
-
-        // The LTE ja LT are chosen such that when A0->A1 is collinear with B0->B1,
-        // or when A1->A2 is collinear with B1->B2, then A is chosen.
-        // This avoids double counting in the case of collinear segments.
-        if (Angle.interiorAngle(a0, intPt, b2) <= bAngle) {
-          area += EdgeVector.area2Term(intPt, a1, a0, false);
-        }
-        if (Angle.interiorAngle(b0, intPt, a2) <= bAngle) {
-          area += EdgeVector.area2Term(intPt, a1, a2, true);
-        }
-        if (Angle.interiorAngle(b0, intPt, a2) < aAngle) {
-          area += EdgeVector.area2Term(intPt, b1, b0, false);
-        }
-        if (Angle.interiorAngle(a0, intPt, b2) < aAngle) {
-          area += EdgeVector.area2Term(intPt, b1, b2, true);
-        }
-      }
-    }
-
-    @Override
-    public boolean isDone() {
-      // Process all intersections
-      return false;
-    }
-  }
-  
 
   private double areaForInteriorVertices(LinearRing ring) {
     /**
diff --git a/modules/lab/src/main/java/org/locationtech/jts/operation/overlayarea/SimpleOverlayArea.java b/modules/lab/src/main/java/org/locationtech/jts/operation/overlayarea/SimpleOverlayArea.java
index 6d386c8af5..00734c91b5 100644
--- a/modules/lab/src/main/java/org/locationtech/jts/operation/overlayarea/SimpleOverlayArea.java
+++ b/modules/lab/src/main/java/org/locationtech/jts/operation/overlayarea/SimpleOverlayArea.java
@@ -11,15 +11,19 @@
  */
 package org.locationtech.jts.operation.overlayarea;
 
-import org.locationtech.jts.algorithm.LineIntersector;
 import org.locationtech.jts.algorithm.Orientation;
 import org.locationtech.jts.algorithm.RayCrossingCounter;
-import org.locationtech.jts.algorithm.RobustLineIntersector;
 import org.locationtech.jts.geom.Coordinate;
 import org.locationtech.jts.geom.CoordinateSequence;
 import org.locationtech.jts.geom.Geometry;
 import org.locationtech.jts.geom.Location;
 import org.locationtech.jts.geom.Polygon;
+import org.locationtech.jts.noding.BasicSegmentString;
+import org.locationtech.jts.noding.SegmentSetMutualIntersector;
+import org.locationtech.jts.noding.SegmentString;
+import org.locationtech.jts.noding.SimpleSegmentSetMutualIntersector;
+
+import java.util.Collections;
 
 /**
  * Computes the result area of an overlay usng the Overlay-Area, for simple polygons.
@@ -78,59 +82,14 @@ private static CoordinateSequence getVertices(Geometry geom) {
   }
   
   private double areaForIntersections(CoordinateSequence ringA, boolean isCCWA, CoordinateSequence ringB, boolean isCCWB) {
-    //TODO: use fast intersection computation?
-    
-    // Compute rays for all intersections
-    LineIntersector li = new RobustLineIntersector();
-    
-    double area = 0;
-    for (int i = 0; i < ringA.size()-1; i++) {
-      Coordinate a0 = ringA.getCoordinate(i);
-      Coordinate a1 = ringA.getCoordinate(i+1);
-      
-      if (isCCWA) {
-        // flip segment orientation
-        Coordinate temp = a0; a0 = a1; a1 = temp;
-      }
-      
-      for (int j = 0; j < ringB.size()-1; j++) {
-        Coordinate b0 = ringB.getCoordinate(j);
-        Coordinate b1 = ringB.getCoordinate(j+1);
-        
-        if (isCCWB) {
-          // flip segment orientation
-          Coordinate temp = b0; b0 = b1; b1 = temp;
-        }
-        
-        li.computeIntersection(a0, a1, b0, b1);
-        if (li.hasIntersection()) {
-          
-          /**
-           * With both rings oriented CW (effectively)
-           * There are two situations for segment intersections:
-           * 
-           * 1) A entering B, B exiting A => rays are IP-A1:R, IP-B0:L
-           * 2) A exiting B, B entering A => rays are IP-A0:L, IP-B1:R
-           * (where :L/R indicates result is to the Left or Right).
-           * 
-           * Use full edge to compute direction, for accuracy.
-           */
-          Coordinate intPt = li.getIntersection(0);
-          
-          boolean isAenteringB = Orientation.COUNTERCLOCKWISE == Orientation.index(a0, a1, b1);
-          
-          if ( isAenteringB ) {
-            area += EdgeVector.area2Term(intPt, a0, a1, true);
-            area += EdgeVector.area2Term(intPt, b1, b0, false);
-          }
-          else {
-            area += EdgeVector.area2Term(intPt, a1, a0, false);
-            area += EdgeVector.area2Term(intPt, b0, b1, true);
-          }
-        }
-      }
-    }
-    return area;
+    SegmentString segStrA = new BasicSegmentString(ringA.toCoordinateArray(), isCCWA);
+    SegmentString segStrB = new BasicSegmentString(ringB.toCoordinateArray(), isCCWB);
+
+    IntersectionVisitor intVisitor = new IntersectionVisitor();
+    SegmentSetMutualIntersector segSetMutInt = new SimpleSegmentSetMutualIntersector(Collections.singleton(segStrA));
+    segSetMutInt.process(Collections.singletonList(segStrB), intVisitor);
+
+    return intVisitor.getArea();
   }
     
   private double areaForInteriorVertices(CoordinateSequence ring, boolean isCCW, CoordinateSequence ring2) {

From bda7a73267c2bdf71244f3e2eba3a450b7ed4d1d Mon Sep 17 00:00:00 2001
From: tanelk <tanel.kiis@gmail.com>
Date: Fri, 22 Mar 2024 22:00:25 +0200
Subject: [PATCH 12/12] Fix SimpleOverlayArea

Signed-off-by: tanelk <tanel.kiis@gmail.com>
---
 .../overlayarea/IntersectionVisitor.java      | 37 ++++++++++++++-----
 1 file changed, 27 insertions(+), 10 deletions(-)

diff --git a/modules/lab/src/main/java/org/locationtech/jts/operation/overlayarea/IntersectionVisitor.java b/modules/lab/src/main/java/org/locationtech/jts/operation/overlayarea/IntersectionVisitor.java
index b81b5b1e4e..4ca5ecfd6f 100644
--- a/modules/lab/src/main/java/org/locationtech/jts/operation/overlayarea/IntersectionVisitor.java
+++ b/modules/lab/src/main/java/org/locationtech/jts/operation/overlayarea/IntersectionVisitor.java
@@ -101,9 +101,14 @@ public void processIntersections(SegmentString a, int aIndex, SegmentString b, i
                 return;
             }
 
-            // If A0->A1 is collinear with B0->B1, then the intersection point from LineIntersector might not be equal to A1 and B1
+            // If A0->A1 is collinear with B0->B1,
+            // then the intersection point from LineIntersector might not be equal to A1 and B1
             Coordinate intPt = a1;
 
+            /* Get the next vertices in the CW direction.
+            Now we have four segments: A0->A1, A1->A2, B0->B1, B1->B2
+            and the intersection point is A1 == B1.
+             */
             Coordinate a2 = a.nextInRing(aIndex + 1);
             Coordinate b2 = b.nextInRing(bIndex + 1);
             if (isCCWA) {
@@ -113,22 +118,34 @@ public void processIntersections(SegmentString a, int aIndex, SegmentString b, i
                 b2 = b.prevInRing(bIndex);
             }
 
-            double aAngle = Angle.interiorAngle(a0, intPt, a2);
-            double bAngle = Angle.interiorAngle(b0, intPt, b2);
+            /* The angles A0->A1->A2 and B0->B1->B2 determine
+             the maximum intersection area interior angle.
+             Edges from the other polygon that lie within this angle
+             are on the boundary of the intersection area.
+
+             Depending on the relative orientation of the polygons,
+             we could pick 0, 2 or 4 segments to contribute to the area.
+
+            The LTE ja LT are chosen such that when A0->A1 is collinear with B0->B1,
+            or when A1->A2 is collinear with B1->B2, then only the segment from polygon A
+            is chosen to avoid double counting.
+             */
+            double angleA0A2 = Angle.interiorAngle(a0, intPt, a2);
+            double angleB0B2 = Angle.interiorAngle(b0, intPt, b2);
+
+            double angleA0B2 = Angle.interiorAngle(a0, intPt, b2);
+            double angleB0A2 = Angle.interiorAngle(b0, intPt, a2);
 
-            // The LTE ja LT are chosen such that when A0->A1 is collinear with B0->B1,
-            // or when A1->A2 is collinear with B1->B2, then A is chosen.
-            // This avoids double counting in the case of collinear segments.
-            if (Angle.interiorAngle(a0, intPt, b2) <= bAngle) {
+            if (angleA0B2 <= angleB0B2) {
                 area += EdgeVector.area2Term(intPt, a1, a0, false);
             }
-            if (Angle.interiorAngle(b0, intPt, a2) <= bAngle) {
+            if (angleB0A2 <= angleB0B2) {
                 area += EdgeVector.area2Term(intPt, a1, a2, true);
             }
-            if (Angle.interiorAngle(b0, intPt, a2) < aAngle) {
+            if (angleB0A2 < angleA0A2) {
                 area += EdgeVector.area2Term(intPt, b1, b0, false);
             }
-            if (Angle.interiorAngle(a0, intPt, b2) < aAngle) {
+            if (angleA0B2 < angleA0A2) {
                 area += EdgeVector.area2Term(intPt, b1, b2, true);
             }
         }