Add MaximumInscribedCircle fast exact calculation for simple shapes (#1123)

Author: dr-jts

modules/core/src/main/java/org/locationtech/jts/algorithm/construct/ExactMaxInscribedCircle.java

@@ -0,0 +1,212 @@
+/*
+ * Copyright (c) 2020 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.algorithm.construct;
+
+import org.locationtech.jts.algorithm.Angle;
+import org.locationtech.jts.algorithm.Orientation;
+import org.locationtech.jts.geom.Coordinate;
+import org.locationtech.jts.geom.CoordinateArrays;
+import org.locationtech.jts.geom.CoordinateSequence;
+import org.locationtech.jts.geom.Geometry;
+import org.locationtech.jts.geom.LineSegment;
+import org.locationtech.jts.geom.LinearRing;
+import org.locationtech.jts.geom.Polygon;
+import org.locationtech.jts.geom.Triangle;
+
+/**
+ * Computes the Maximum Inscribed Circle for some kinds of convex polygons.
+ * It determines the circle center point by computing Voronoi node points
+ * and testing them for distance to generating edges.
+ * This is more precise than iterated approximation, 
+ * and faster for small polygons (such as triangles and convex quadrilaterals).
+ * 
+ * @author Martin Davis
+ *
+ */
+class ExactMaxInscribedCircle {
+  
+  /**
+   * Tests whether a given geometry is supported by this class.
+   * Currently only triangles and convex quadrilaterals are supported.
+   * 
+   * @param geom an areal geometry
+   * @return true if the geometry shape can be evaluated
+   */
+  public static boolean isSupported(Geometry geom) {
+    if (! isSimplePolygon(geom)) 
+      return false;
+    Polygon polygon = (Polygon) geom;
+    if (isTriangle(polygon)) 
+      return true;
+    if (isQuadrilateral(polygon) && isConvex(polygon))
+      return true;
+    return false;
+  }
+
+  private static boolean isSimplePolygon(Geometry geom) {
+    return geom instanceof Polygon
+        && ((Polygon) geom).getNumInteriorRing() == 0; 
+  }
+
+  private static boolean isTriangle(Polygon polygon) {
+    return polygon.getNumPoints() == 4;
+  }
+  
+  private static boolean isQuadrilateral(Polygon polygon) {
+    return polygon.getNumPoints() == 5;
+  }
+
+  public static Coordinate[] computeRadius(Polygon polygon) {
+    Coordinate[] ring = polygon.getExteriorRing().getCoordinates();
+    if (ring.length == 4)
+      return computeTriangle(ring);
+    else if (ring.length == 5)
+      return computeConvexQuadrilateral(ring);
+    throw new IllegalArgumentException("Input must be a triangle or convex quadrilateral");
+  }
+  
+  private static Coordinate[] computeTriangle(Coordinate[] ring) {
+    Coordinate center = Triangle.inCentre(ring[0], ring[1], ring[2]);
+    LineSegment seg = new LineSegment(ring[0], ring[1]);
+    Coordinate radius = seg.project(center);
+    return new Coordinate[] { center, radius };
+  }
+
+  /**
+   * The Voronoi nodes of a convex polygon occur at the intersection point
+   * of two bisectors of each triplet of edges.
+   * The Maximum Inscribed Circle center is the node
+   * is the farthest distance from the generating edges.
+   * For a quadrilateral there are 4 distinct edge triplets, 
+   * at each edge with its adjacent edges. 
+   * 
+   * @param ring the polygon ring
+   * @return an array containing the incircle center and radius points
+   */
+  private static Coordinate[] computeConvexQuadrilateral(Coordinate[] ring) {
+    Coordinate[] ringCW = CoordinateArrays.orient(ring, true);
+    
+    double diameter = CoordinateArrays.envelope(ringCW).getDiameter();
+    
+    //-- compute corner bisectors
+    LineSegment[] bisector = computeBisectors(ringCW, diameter);
+    //-- compute nodes and find interior one farthest from sides
+    double maxDist = -1;
+    Coordinate center = null;
+    Coordinate radius = null;
+    for (int i = 0; i < 4; i++) {
+      LineSegment b1 = bisector[i];
+      int i2 = (i + 1) % 4;
+      LineSegment b2 = bisector[i2];
+
+      Coordinate nodePt = b1.intersection(b2);
+      
+      //-- only interior nodes are considered
+      if (! isPointInConvexRing(ringCW, nodePt)) {
+        continue;
+      }
+      
+      //-- check if node is further than current max center
+      Coordinate r = nearestEdgePt(ringCW, nodePt);
+      double dist = nodePt.distance(r);
+      if (maxDist < 0 || dist > maxDist) {
+        center = nodePt;
+        radius = r;
+        maxDist = dist;
+        //System.out.println(WKTWriter.toLineString(center, radius));
+      }
+    }
+    return new Coordinate[] { center, radius };
+  }
+
+  private static LineSegment[] computeBisectors(Coordinate[] ptsCW, double diameter) {
+    LineSegment[] bisector = new LineSegment[4];
+    for (int i = 0; i < 4; i++) {
+      bisector[i] = computeConvexBisector(ptsCW, i, diameter);
+    }
+    return bisector;
+  }
+
+  private static Coordinate nearestEdgePt(Coordinate[] ring, Coordinate pt) {
+    Coordinate nearestPt = null;
+    double minDist = -1;
+    for (int i = 0; i < ring.length - 1; i++) {
+      LineSegment edge = new LineSegment(ring[i], ring[i + 1]);
+      Coordinate r = edge.closestPoint(pt);
+      double dist = pt.distance(r);
+      if (minDist < 0 || dist < minDist) {
+        minDist = dist;
+        nearestPt = r;
+      }
+    }
+    return nearestPt;
+  }
+
+  private static LineSegment computeConvexBisector(Coordinate[] pts, int index, double len) {
+    Coordinate basePt = pts[index];
+    int iPrev = index == 0 ? pts.length - 2 : index - 1;
+    int iNext = index >= pts.length ? 0 : index + 1;
+    Coordinate pPrev = pts[iPrev];
+    Coordinate pNext = pts[iNext];
+    if (! isConvex(pPrev, basePt, pNext))
+      throw new IllegalArgumentException("Input is not convex");
+    double bisectAng = Angle.bisector(pPrev, basePt, pNext);
+    Coordinate endPt = Angle.project(basePt, bisectAng, len);
+    return new LineSegment(basePt.copy(), endPt);
+  }
+  
+  private static boolean isConvex(Polygon polygon) {
+    LinearRing shell = polygon.getExteriorRing();
+    return isConvex(shell.getCoordinateSequence());
+  }
+  
+  private static boolean isConvex(CoordinateSequence ring) {
+    /**
+     * A ring cannot be all concave, so if it has a consistent
+     * orientation it must be convex.
+     */
+    int n = ring.size();
+    if (n < 4) 
+      return false;
+    int ringOrient = 0;
+    for (int i = 0; i < n - 1; i++) {
+      int i1 = i + 1;
+      int i2 = (i1 >= n - 1) ? 1 : i1 + 1;
+      int orient = Orientation.index(ring.getCoordinate(i), 
+          ring.getCoordinate(i1), ring.getCoordinate(i2));
+      if (orient == Orientation.COLLINEAR)
+        continue;
+      if (ringOrient == 0) {
+        ringOrient = orient;
+      }
+      else if (orient != ringOrient) {
+          return false;
+      }
+    }
+    return true;
+  }
+
+  private static boolean isConvex(Coordinate p0, Coordinate p1, Coordinate p2) {
+    return Orientation.CLOCKWISE == Orientation.index(p0, p1, p2);
+  }
+
+  private static boolean isPointInConvexRing(Coordinate[] ringCW, Coordinate p) {
+    for (int i = 0; i < ringCW.length - 1; i++) {
+      Coordinate p0 = ringCW[i];
+      Coordinate p1 = ringCW[i + 1];
+      int orient = Orientation.index(p0, p1, p);
+      if (orient == Orientation.COUNTERCLOCKWISE)
+        return false;
+    }
+    return true;
+  }
+}

modules/core/src/main/java/org/locationtech/jts/algorithm/construct/MaximumInscribedCircle.java

@@ -144,8 +144,6 @@ public MaximumInscribedCircle(Geometry polygonal, double tolerance) {
     this.inputGeom = polygonal;
     this.factory = polygonal.getFactory();
     this.tolerance = tolerance;
-    ptLocater = new IndexedPointInAreaLocator(polygonal);
-    indexedDistance = new IndexedFacetDistance( polygonal.getBoundary() );
   }
 
   /**
@@ -210,7 +208,39 @@ private double distanceToBoundary(double x, double y) {
 
   private void compute() {
     // check if already computed
-    if (centerCell != null) return;
+    if (centerPt != null) return;
+    
+    /**
+     * Handle empty or flat geometries.
+     */
+    if (inputGeom.getArea() == 0.0) {
+      Coordinate c = inputGeom.getCoordinate().copy();
+      createResult(c, c.copy());
+      return;
+    }
+    
+    /**
+     * Optimization for small simple convex polygons 
+     */
+    if (ExactMaxInscribedCircle.isSupported(inputGeom)) {
+      Coordinate[] centreRadius = ExactMaxInscribedCircle.computeRadius((Polygon) inputGeom);
+      createResult(centreRadius[0], centreRadius[1]);
+      return;
+    }
+    
+    computeApproximation();
+  }
+
+  private void createResult(Coordinate c, Coordinate r) {
+    centerPt = c;
+    radiusPt = r;
+    centerPoint = factory.createPoint(centerPt);
+    radiusPoint = factory.createPoint(radiusPt);
+  }
+
+  private void computeApproximation() {  
+    ptLocater = new IndexedPointInAreaLocator(inputGeom);
+    indexedDistance = new IndexedFacetDistance( inputGeom.getBoundary() );
 
     // Priority queue of cells, ordered by maximum distance from boundary
     PriorityQueue<Cell> cellQueue = new PriorityQueue<>();

modules/core/src/main/java/org/locationtech/jts/geom/CoordinateArrays.java

@@ -15,6 +15,7 @@
 import java.util.Collection;
 import java.util.Comparator;
 
+import org.locationtech.jts.algorithm.Orientation;
 import org.locationtech.jts.math.MathUtil;
 
 
@@ -192,6 +193,24 @@ public static Coordinate ptNotInList(Coordinate[] testPts, Coordinate[] pts) {
     return null;
   }
 
+  /**
+   * Orients an array of points to have a specified orientation.
+   * The array is not copied if it already has the required orientation.
+   * The points are not copied.
+   * 
+   * @param pts the array to orient
+   * @param orientCW true if the points should be oriented CVW
+   * @return the oriented array
+   */
+  public static Coordinate[] orient(Coordinate[] pts, boolean orientCW) {
+    boolean isFlipped = orientCW == Orientation.isCCW(pts);
+    if (isFlipped) {
+      pts = pts.clone();
+      CoordinateArrays.reverse(pts);
+    }
+    return pts;
+  }
+  
   /**
    * Compares two {@link Coordinate} arrays
    * in the forward direction of their coordinates,

modules/core/src/main/java/org/locationtech/jts/geom/Triangle.java

@@ -261,7 +261,8 @@ private static double det(double m00, double m01, double m10, double m11)
    * the point which is equidistant from the sides of the triangle. It is also
    * the point at which the bisectors of the triangle's angles meet. It is the
    * centre of the triangle's <i>incircle</i>, which is the unique circle that
-   * is tangent to each of the triangle's three sides.
+   * is tangent to each of the triangle's three sides
+   * (and hence the maximum inscribed circle).
    * <p>
    * The incentre always lies within the triangle.
    * 

modules/core/src/test/java/org/locationtech/jts/algorithm/construct/MaximumInscribedCircleTest.java

@@ -15,16 +15,41 @@ public static void main(String args[]) {
 
   public MaximumInscribedCircleTest(String name) { super(name); }
 
+  public void testTriangleRight() {
+    checkCircle("POLYGON ((1 1, 1 9, 9 1, 1 1))", 
+       0.001, 3.343, 3.343, 2.343 );
+  }
+
+  public void testTriangleObtuse() {
+    checkCircle("POLYGON ((1 1, 1 9, 2 2, 1 1))", 
+       0.001, 1.485, 2.173, 0.485 );
+  }
+  
   public void testSquare() {
     checkCircle("POLYGON ((100 200, 200 200, 200 100, 100 100, 100 200))", 
        0.001, 150, 150, 50 );
   }
 
+  public void testThinQuad() {
+    checkCircle("POLYGON ((1 2, 9 3, 9 1, 1 1, 1 2))", 
+       0.001, 8.0623, 1.9377, 0.93774 );
+  }
+
   public void testDiamond() {
     checkCircle("POLYGON ((150 250, 50 150, 150 50, 250 150, 150 250))", 
        0.001, 150, 150, 70.71 );
   }
 
+  public void testChevron() {
+    checkCircle("POLYGON ((1 1, 6 9, 3.7 2.5, 9 1, 1 1))", 
+       0.001, 2.82, 2.008, 1.008 );
+  }
+
+  public void testChevronFat() {
+    checkCircle("POLYGON ((1 1, 6 9, 5.9 5, 9 1, 1 1))", 
+       0.001, 4.7545, 3.0809, 2.081 );
+  }
+
   public void testCircle() {
     Geometry centre = read("POINT (100 100)");
     Geometry circle = centre.buffer(100, 20);
@@ -65,7 +90,7 @@ public void testCollapsedLineFlat() {
   }
 
   /**
-   * Invalid polygon collapsed to a point
+   * Invalid triangle polygon collapsed to a point
    */
   public void testCollapsedPoint() {
     checkCircle("POLYGON ((100 100, 100 100, 100 100, 100 100))", 

modules/core/src/test/java/org/locationtech/jts/geom/CoordinateArraysTest.java

@@ -12,6 +12,8 @@
 package org.locationtech.jts.geom;
 
 
+import org.locationtech.jts.algorithm.Orientation;
+
 import junit.textui.TestRunner;
 import test.jts.GeometryTestCase;
 
@@ -174,6 +176,31 @@ public void testEnforceConsistency(){
     assertTrue( array[1] != fixed[1] ); // processing needed to CoordinateXYZM
   }
 
+  public void testOrientCW() {
+    checkOrient("POLYGON ((1 1, 9 9, 9 1, 1 1))");
+  }
+  
+  public void testOrientCCW() {
+    checkOrient("POLYGON ((9 7, 5 9, 1 4, 5 4, 4 1, 8 1, 9 7))");
+  }
+  
+  private void checkOrient(String wkt) {
+    Coordinate[] pts = read(wkt).getCoordinates();
+    //-- orient CW
+    Coordinate[] ptsCW = CoordinateArrays.orient(pts, true);
+    assertEquals(false, Orientation.isCCW(ptsCW));
+    Coordinate[] ptsCCW = CoordinateArrays.orient(pts, false);
+    assertEquals(true, Orientation.isCCW(ptsCCW));
+    //-- check that original is unchanged for same orientation
+    boolean isCCW = Orientation.isCCW(pts);
+    if (isCCW) {
+      assertTrue(pts == ptsCCW);
+    }
+    else {
+      assertTrue(pts == ptsCW);      
+    }
+  }
+
   private static void checkCoordinateAt(Coordinate[] seq1, int pos1,
                                         Coordinate[] seq2, int pos2) {
     Coordinate c1 = seq1[pos1], c2 = seq2[pos2];