Fix vertex intersection in OverlayArea

modules/lab/src/main/java/org/locationtech/jts/operation/overlayarea/IntersectionVisitor.java

@@ -0,0 +1,159 @@
+/*
+ * 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;
+
+            /* 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) {
+                a2 = a.prevInRing(aIndex);
+            }
+            if (isCCWB) {
+                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 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);
+
+            if (angleA0B2 <= angleB0B2) {
+                area += EdgeVector.area2Term(intPt, a1, a0, false);
+            }
+            if (angleB0A2 <= angleB0B2) {
+                area += EdgeVector.area2Term(intPt, a1, a2, true);
+            }
+            if (angleB0A2 < angleA0A2) {
+                area += EdgeVector.area2Term(intPt, b1, b0, false);
+            }
+            if (angleA0B2 < angleA0A2) {
+                area += EdgeVector.area2Term(intPt, b1, b2, true);
+            }
+        }
+    }
+
+    @Override
+    public boolean isDone() {
+        // Process all intersections
+        return false;
+    }
+}

modules/lab/src/main/java/org/locationtech/jts/operation/overlayarea/OverlayArea.java

@@ -11,12 +11,8 @@
  */
 package org.locationtech.jts.operation.overlayarea;
 
-import java.util.List;
-
 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;
@@ -25,21 +21,25 @@
 import org.locationtech.jts.geom.Envelope;
 import org.locationtech.jts.geom.Geometry;
 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.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.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
  * 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.
@@ -74,12 +74,10 @@ 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;
-  private STRtree indexSegs;
+  private SegmentSetMutualIntersector segSetMutInt;
   private KdTree vertexIndex;
 
   public OverlayArea(Geometry geom) {
@@ -92,7 +90,7 @@ public OverlayArea(Geometry geom) {
 
     geomEnv0 = geom.getEnvelopeInternal();
     locator0 = new IndexedPointInAreaLocator(geom);
-    indexSegs = buildSegmentIndex(geom);
+    segSetMutInt = buildSegmentIndex(geom);
     vertexIndex = buildVertexIndex(geom);
   }
 
@@ -173,7 +171,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;
@@ -207,80 +205,16 @@ 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));
 
-  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;
-    }
-
-    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);
-    }
-    else {
-      return EdgeVector.area2Term(intPt, a1, a0, false)
-       + EdgeVector.area2Term(intPt, b0, b1, true);
-    }
+    IntersectionVisitor intVisitor = new IntersectionVisitor();
+    segSetMutInt.process(Collections.singletonList(segStr), intVisitor);
+
+    return intVisitor.getArea();
   }
-
+
+
   private double areaForInteriorVertices(LinearRing ring) {
     /**
      * Compute rays originating at vertices inside the intersection result
@@ -335,22 +269,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) {