Skip to content

Commit

Permalink
Add loop line merger (#1083)
Browse files Browse the repository at this point in the history
  • Loading branch information
wipfli authored Nov 25, 2024
1 parent 8a1b54f commit c470dc3
Show file tree
Hide file tree
Showing 11 changed files with 1,338 additions and 22 deletions.
Original file line number Diff line number Diff line change
@@ -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));
}
}
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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());
Expand All @@ -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);
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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;
Expand Down
Loading

0 comments on commit c470dc3

Please sign in to comment.