Add loop line merger

planetiler-benchmarks/src/main/java/com/onthegomap/planetiler/benchmarks/BenchmarkLineMerge.java

@@ -0,0 +1,114 @@
+package com.onthegomap.planetiler.benchmarks;
+
+import com.onthegomap.planetiler.geo.GeoUtils;
+import com.onthegomap.planetiler.util.Format;
+import com.onthegomap.planetiler.util.FunctionThatThrows;
+import com.onthegomap.planetiler.util.Gzip;
+import com.onthegomap.planetiler.util.LoopLineMerger;
+import java.io.IOException;
+import java.math.BigDecimal;
+import java.math.MathContext;
+import java.nio.file.Files;
+import java.nio.file.Path;
+import java.time.Duration;
+import java.util.ArrayList;
+import java.util.Collection;
+import java.util.List;
+import org.locationtech.jts.geom.Coordinate;
+import org.locationtech.jts.geom.CoordinateXY;
+import org.locationtech.jts.geom.Geometry;
+import org.locationtech.jts.geom.GeometryFactory;
+import org.locationtech.jts.geom.LineString;
+import org.locationtech.jts.io.ParseException;
+import org.locationtech.jts.io.WKBReader;
+import org.locationtech.jts.operation.linemerge.LineMerger;
+
+public class BenchmarkLineMerge {
+  private static int numLines;
+
+  public static void main(String[] args) throws Exception {
+    for (int i = 0; i < 10; i++) {
+      time("       JTS", geom -> {
+        var lm = new LineMerger();
+        lm.add(geom);
+        return lm.getMergedLineStrings();
+      });
+      time("   loop(0)", geom -> loopMerger(0).add(geom).getMergedLineStrings());
+      time(" loop(0.1)", geom -> loopMerger(0.1).add(geom).getMergedLineStrings());
+      time("loop(20.0)", geom -> loopMerger(20).add(geom).getMergedLineStrings());
+    }
+    System.err.println(numLines);
+  }
+
+  private static LoopLineMerger loopMerger(double minLength) {
+    var lm = new LoopLineMerger();
+    lm.setMinLength(minLength);
+    lm.setStubMinLength(minLength);
+    lm.setLoopMinLength(minLength);
+    lm.setTolerance(1);
+    lm.setMergeStrokes(true);
+    return lm;
+  }
+
+  private static void time(String name, FunctionThatThrows<Geometry, Collection<LineString>> fn) throws Exception {
+    System.err.println(String.join("\t",
+      name,
+      timeMillis(read("mergelines_200433_lines.wkb.gz"), fn),
+      timeMillis(read("mergelines_239823_lines.wkb.gz"), fn),
+      "(/s):",
+      timePerSec(read("mergelines_1759_point_line.wkb.gz"), fn),
+      timePerSec(makeLines(50, 2), fn),
+      timePerSec(makeLines(10, 10), fn),
+      timePerSec(makeLines(2, 50), fn)
+    ));
+  }
+
+  private static String timePerSec(Geometry geometry, FunctionThatThrows<Geometry, Collection<LineString>> fn)
+    throws Exception {
+    long start = System.nanoTime();
+    long end = start + Duration.ofSeconds(1).toNanos();
+    int num = 0;
+    for (; System.nanoTime() < end;) {
+      numLines += fn.apply(geometry).size();
+      num++;
+    }
+    return Format.defaultInstance()
+      .numeric(Math.round(num * 1d / ((System.nanoTime() - start) * 1d / Duration.ofSeconds(1).toNanos())), true);
+  }
+
+  private static String timeMillis(Geometry geometry, FunctionThatThrows<Geometry, Collection<LineString>> fn)
+    throws Exception {
+    long start = System.nanoTime();
+    long end = start + Duration.ofSeconds(1).toNanos();
+    int num = 0;
+    for (; System.nanoTime() < end;) {
+      numLines += fn.apply(geometry).size();
+      num++;
+    }
+    // equivalent of toPrecision(3)
+    long nanosPer = (System.nanoTime() - start) / num;
+    var bd = new BigDecimal(nanosPer, new MathContext(3));
+    return Format.padRight(Duration.ofNanos(bd.longValue()).toString().replace("PT", ""), 6);
+  }
+
+
+  private static Geometry read(String fileName) throws IOException, ParseException {
+    var path = Path.of("planetiler-core", "src", "test", "resources", "mergelines", fileName);
+    byte[] bytes = Gzip.gunzip(Files.readAllBytes(path));
+    return new WKBReader().read(bytes);
+  }
+
+  private static Geometry makeLines(int lines, int parts) {
+    List<LineString> result = new ArrayList<>();
+    double idx = 0;
+    for (int i = 0; i < lines; i++) {
+      Coordinate[] coords = new Coordinate[parts];
+      for (int j = 0; j < parts; j++) {
+        coords[j] = new CoordinateXY(idx, idx);
+        idx += 0.5;
+      }
+      result.add(GeoUtils.JTS_FACTORY.createLineString(coords));
+    }
+    return new GeometryFactory().createMultiLineString(result.toArray(LineString[]::new));
+  }
+}

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

@@ -4,13 +4,13 @@
 import com.carrotsearch.hppc.IntObjectMap;
 import com.carrotsearch.hppc.IntStack;
 import com.onthegomap.planetiler.collection.Hppc;
-import com.onthegomap.planetiler.geo.DouglasPeuckerSimplifier;
 import com.onthegomap.planetiler.geo.GeoUtils;
 import com.onthegomap.planetiler.geo.GeometryException;
 import com.onthegomap.planetiler.geo.GeometryType;
 import com.onthegomap.planetiler.geo.MutableCoordinateSequence;
 import com.onthegomap.planetiler.stats.DefaultStats;
 import com.onthegomap.planetiler.stats.Stats;
+import com.onthegomap.planetiler.util.LoopLineMerger;
 import java.util.ArrayList;
 import java.util.BitSet;
 import java.util.Collection;
@@ -171,7 +171,12 @@ public static List<VectorTile.Feature> mergeLineStrings(List<VectorTile.Feature>
       if (groupedFeatures.size() == 1 && buffer == 0d && lengthLimit == 0 && (!resimplify || tolerance == 0)) {
         result.add(feature1);
       } else {
-        LineMerger merger = new LineMerger();
+        LoopLineMerger merger = new LoopLineMerger()
+          .setTolerance(tolerance)
+          .setMergeStrokes(true)
+          .setMinLength(lengthLimit)
+          .setLoopMinLength(lengthLimit)
+          .setStubMinLength(0.5);
         for (VectorTile.Feature feature : groupedFeatures) {
           try {
             merger.add(feature.geometry().decode());
@@ -180,24 +185,14 @@ public static List<VectorTile.Feature> mergeLineStrings(List<VectorTile.Feature>
           }
         }
         List<LineString> outputSegments = new ArrayList<>();
-        for (Object merged : merger.getMergedLineStrings()) {
-          if (merged instanceof LineString line && line.getLength() >= lengthLimit) {
-            // re-simplify since some endpoints of merged segments may be unnecessary
-            if (line.getNumPoints() > 2 && tolerance >= 0) {
-              Geometry simplified = DouglasPeuckerSimplifier.simplify(line, tolerance);
-              if (simplified instanceof LineString simpleLineString) {
-                line = simpleLineString;
-              } else {
-                LOGGER.warn("line string merge simplify emitted {}", simplified.getGeometryType());
-              }
-            }
-            if (buffer >= 0) {
-              removeDetailOutsideTile(line, buffer, outputSegments);
-            } else {
-              outputSegments.add(line);
-            }
+        for (var line : merger.getMergedLineStrings()) {
+          if (buffer >= 0) {
+            removeDetailOutsideTile(line, buffer, outputSegments);
+          } else {
+            outputSegments.add(line);
           }
         }
+
         if (!outputSegments.isEmpty()) {
           outputSegments = sortByHilbertIndex(outputSegments);
           Geometry newGeometry = GeoUtils.combineLineStrings(outputSegments);

planetiler-core/src/main/java/com/onthegomap/planetiler/geo/DouglasPeuckerSimplifier.java

@@ -11,6 +11,8 @@
  */
 package com.onthegomap.planetiler.geo;
 
+import java.util.ArrayList;
+import java.util.List;
 import org.locationtech.jts.geom.Coordinate;
 import org.locationtech.jts.geom.CoordinateSequence;
 import org.locationtech.jts.geom.Geometry;
@@ -45,6 +47,22 @@ public static Geometry simplify(Geometry geom, double distanceTolerance) {
     return (new DPTransformer(distanceTolerance)).transform(geom);
   }
 
+  /**
+   * Returns a copy of {@code coords}, simplified using Douglas Peucker Algorithm.
+   *
+   * @param coords            the coordinate list to simplify
+   * @param distanceTolerance the threshold below which we discard points
+   * @param area              true if this is a polygon to retain at least 4 points to avoid collapse
+   * @return the simplified coordinate list
+   */
+  public static List<Coordinate> simplify(List<Coordinate> coords, double distanceTolerance, boolean area) {
+    if (coords.isEmpty()) {
+      return List.of();
+    }
+
+    return (new DPTransformer(distanceTolerance)).transformCoordinateList(coords, area);
+  }
+
   private static class DPTransformer extends GeometryTransformer {
 
     private final double sqTolerance;
@@ -84,6 +102,42 @@ private static double getSqSegDist(double px, double py, double p1x, double p1y,
       return dx * dx + dy * dy;
     }
 
+    private void subsimplify(List<Coordinate> in, List<Coordinate> out, int first, int last, int numForcedPoints) {
+      // numForcePoints lets us keep some points even if they are below simplification threshold
+      boolean force = numForcedPoints > 0;
+      double maxSqDist = force ? -1 : sqTolerance;
+      int index = -1;
+      Coordinate p1 = in.get(first);
+      Coordinate p2 = in.get(last);
+      double p1x = p1.x;
+      double p1y = p1.y;
+      double p2x = p2.x;
+      double p2y = p2.y;
+
+      int i = first + 1;
+      Coordinate furthest = null;
+      for (Coordinate coord : in.subList(first + 1, last)) {
+        double sqDist = getSqSegDist(coord.x, coord.y, p1x, p1y, p2x, p2y);
+
+        if (sqDist > maxSqDist) {
+          index = i;
+          furthest = coord;
+          maxSqDist = sqDist;
+        }
+        i++;
+      }
+
+      if (force || maxSqDist > sqTolerance) {
+        if (index - first > 1) {
+          subsimplify(in, out, first, index, numForcedPoints - 1);
+        }
+        out.add(furthest);
+        if (last - index > 1) {
+          subsimplify(in, out, index, last, numForcedPoints - 2);
+        }
+      }
+    }
+
     private void subsimplify(CoordinateSequence in, MutableCoordinateSequence out, int first, int last,
       int numForcedPoints) {
       // numForcePoints lets us keep some points even if they are below simplification threshold
@@ -117,6 +171,20 @@ private void subsimplify(CoordinateSequence in, MutableCoordinateSequence out, i
       }
     }
 
+    protected List<Coordinate> transformCoordinateList(List<Coordinate> coords, boolean area) {
+      if (coords.isEmpty()) {
+        return coords;
+      }
+      // make sure we include the first and last points even if they are closer than the simplification threshold
+      List<Coordinate> result = new ArrayList<>();
+      result.add(coords.getFirst());
+      // for polygons, additionally keep at least 2 intermediate points even if they are below simplification threshold
+      // to avoid collapse.
+      subsimplify(coords, result, 0, coords.size() - 1, area ? 2 : 0);
+      result.add(coords.getLast());
+      return result;
+    }
+
     @Override
     protected CoordinateSequence transformCoordinates(CoordinateSequence coords, Geometry parent) {
       boolean area = parent instanceof LinearRing;