Add API for maximum inscribed circle/pole of inaccessibility centerpoint of a polygon (#723)

planetiler-core/src/main/java/com/onthegomap/planetiler/FeatureCollector.java

@@ -15,6 +15,7 @@ import java.util.Iterator;
 import java.util.List;
 import java.util.Map;
 import java.util.TreeMap;
+import org.locationtech.jts.algorithm.construct.MaximumInscribedCircle;
 import org.locationtech.jts.geom.Geometry;
 
 /**
@@ -172,6 +173,33 @@ public class FeatureCollector implements Iterable<FeatureCollector.Feature> {
     }
   }
 
+  /**
+   * Starts building a new point map feature at the furthest interior point of a polygon from its edge using
+   * {@link MaximumInscribedCircle} (aka "pole of inaccessibility") of the source feature.
+   * <p>
+   * NOTE: This is substantially more expensive to compute than {@link #centroid(String)} or
+   * {@link #pointOnSurface(String)}, especially for small {@code tolerance} values.
+   *
+   * @param layer     the output vector tile layer this feature will be written to
+   * @param tolerance precision for calculating maximum inscribed circle. 0.01 means 1% of the square root of the area.
+   *                  Smaller values for a more precise tolerance become very expensive to compute. Values between 5%
+   *                  and 10% are a good compromise of performance vs. precision.
+   * @return a feature that can be configured further.
+   */
+  public Feature innermostPoint(String layer, double tolerance) {
+    try {
+      return geometry(layer, source.innermostPoint(tolerance));
+    } catch (GeometryException e) {
+      e.log(stats, "feature_innermost_point", "Error constructing innermost point for " + source.id());
+      return new Feature(layer, EMPTY_GEOM, source.id());
+    }
+  }
+
+  /** Alias for {@link #innermostPoint(String, double)} with a default tolerance of 10%. */
+  public Feature innermostPoint(String layer) {
+    return innermostPoint(layer, 0.1);
+  }
+
   /**
    * Creates new feature collector instances for each source feature that we encounter.
    */
@@ -263,8 +291,8 @@ public class FeatureCollector implements Iterable<FeatureCollector.Feature> {
 
     /**
      * Sets the value by which features are sorted within a layer in the output vector tile. Sort key gets packed into
-     * {@link FeatureGroup#SORT_KEY_BITS} bits so the range of this is limited to {@code -(2^(bits-1))} to {@code
-     * (2^(bits-1))-1}.
+     * {@link FeatureGroup#SORT_KEY_BITS} bits so the range of this is limited to {@code -(2^(bits-1))} to
+     * {@code (2^(bits-1))-1}.
      * <p>
      * Circles, lines, and polygons are rendered in the order they appear in each layer, so features that appear later
      * (higher sort key) show up on top of features with a lower sort key.

planetiler-core/src/main/java/com/onthegomap/planetiler/reader/SourceFeature.java

@@ -7,6 +7,7 @@ import com.onthegomap.planetiler.reader.osm.OsmRelationInfo;
 import java.util.ArrayList;
 import java.util.List;
 import java.util.Map;
+import org.locationtech.jts.algorithm.construct.MaximumInscribedCircle;
 import org.locationtech.jts.geom.Geometry;
 import org.locationtech.jts.geom.LineString;
 import org.locationtech.jts.geom.Lineal;
@@ -34,6 +35,8 @@ public abstract class SourceFeature implements WithTags, WithGeometryType {
   private Geometry centroid = null;
   private Geometry pointOnSurface = null;
   private Geometry centroidIfConvex = null;
+  private double innermostPointTolerance = Double.NaN;
+  private Geometry innermostPoint = null;
   private Geometry linearGeometry = null;
   private Geometry polygonGeometry = null;
   private Geometry validPolygon = null;
@@ -126,6 +129,26 @@ public abstract class SourceFeature implements WithTags, WithGeometryType {
         worldGeometry().getInteriorPoint());
   }
 
+  /**
+   * Returns {@link MaximumInscribedCircle#getCenter()} of this geometry in world web mercator coordinates.
+   *
+   * @param tolerance precision for calculating maximum inscribed circle. 0.01 means 1% of the square root of the area.
+   *                  Smaller values for a more precise tolerance become very expensive to compute. Values between
+   *                  0.05-0.1 are a good compromise of performance vs. precision.
+   */
+  public final Geometry innermostPoint(double tolerance) throws GeometryException {
+    if (canBePolygon()) {
+      // cache as long as the tolerance hasn't changed
+      if (tolerance != innermostPointTolerance || innermostPoint == null) {
+        innermostPoint = MaximumInscribedCircle.getCenter(polygon(), Math.sqrt(area()) * tolerance);
+        innermostPointTolerance = tolerance;
+      }
+      return innermostPoint;
+    } else {
+      return pointOnSurface();
+    }
+  }
+
   private Geometry computeCentroidIfConvex() throws GeometryException {
     if (!canBePolygon()) {
       return centroid();

planetiler-core/src/test/java/com/onthegomap/planetiler/FeatureCollectorTest.java

@@ -493,6 +493,34 @@ class FeatureCollectorTest {
     assertFalse(iter.hasNext());
   }
 
+  @Test
+  void testInnermostPoint() {
+    /*
+      _____
+     | · __|
+     |__|
+     */
+    var sourceLine = newReaderFeature(newPolygon(worldToLatLon(
+      0, 0,
+      1, 0,
+      1, 0.5,
+      0.5, 0.5,
+      0.5, 1,
+      0, 1,
+      0, 0
+    )), Map.of());
+
+    var fc = factory.get(sourceLine);
+    fc.innermostPoint("layer").setZoomRange(0, 10);
+    var iter = fc.iterator();
+
+    var item = iter.next();
+    assertEquals(GeometryType.POINT, item.getGeometryType());
+    assertEquals(round(newPoint(0.28, 0.28)), round(item.getGeometry(), 1e2));
+
+    assertFalse(iter.hasNext());
+  }
+
   @Test
   void testMultiPolygonCoercion() throws GeometryException {
     var sourceLine = newReaderFeature(newMultiPolygon(