progress and minheap proto

pull/1109/head
Mike Barry 2024-11-21 06:12:42 -05:00
rodzic 00b0c8bec6
commit bf3dcb3d85
5 zmienionych plików z 803 dodań i 0 usunięć

Wyświetl plik

@ -0,0 +1,84 @@
package com.onthegomap.planetiler.benchmarks;
import com.onthegomap.planetiler.geo.DouglasPeuckerSimplifier;
import com.onthegomap.planetiler.geo.VWSimplifier;
import com.onthegomap.planetiler.util.Format;
import com.onthegomap.planetiler.util.FunctionThatThrows;
import java.math.BigDecimal;
import java.math.MathContext;
import java.time.Duration;
import org.locationtech.jts.geom.CoordinateXY;
import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.util.GeometricShapeFactory;
public class BenchmarkSimplify {
private static int numLines;
public static void main(String[] args) throws Exception {
for (int i = 0; i < 10; i++) {
time(" DP(0.1)", geom -> DouglasPeuckerSimplifier.simplify(geom, 0.1));
time(" DP(1)", geom -> DouglasPeuckerSimplifier.simplify(geom, 1));
time(" DP(20)", geom -> DouglasPeuckerSimplifier.simplify(geom, 20));
time(" JTS VW(0)", geom -> org.locationtech.jts.simplify.VWSimplifier.simplify(geom, 0.01));
time("JTS VW(0.1)", geom -> org.locationtech.jts.simplify.VWSimplifier.simplify(geom, 0.1));
time(" JTS VW(1)", geom -> org.locationtech.jts.simplify.VWSimplifier.simplify(geom, 1));
time(" JTS VW(20)", geom -> org.locationtech.jts.simplify.VWSimplifier.simplify(geom, 20));
time(" VW(0)", geom -> new VWSimplifier().setTolerance(0).setWeight(0.7).transform(geom));
time(" VW(0.1)", geom -> new VWSimplifier().setTolerance(0.1).setWeight(0.7).transform(geom));
time(" VW(1)", geom -> new VWSimplifier().setTolerance(1).setWeight(0.7).transform(geom));
time(" VW(20)", geom -> new VWSimplifier().setTolerance(20).setWeight(0.7).transform(geom));
}
System.err.println(numLines);
}
private static void time(String name, FunctionThatThrows<Geometry, Geometry> fn) throws Exception {
System.err.println(String.join("\t",
name,
timePerSec(makeLines(2), fn),
timePerSec(makeLines(10), fn),
timePerSec(makeLines(50), fn),
timePerSec(makeLines(100), fn),
timePerSec(makeLines(10_000), fn)
));
}
private static String timePerSec(Geometry geometry, FunctionThatThrows<Geometry, Geometry> fn)
throws Exception {
long start = System.nanoTime();
long end = start + Duration.ofSeconds(1).toNanos();
int num = 0;
boolean first = true;
for (; System.nanoTime() < end;) {
numLines += fn.apply(geometry).getNumPoints();
if (first) {
first = false;
}
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, Geometry> 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).getNumPoints();
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 makeLines(int parts) {
var shapeFactory = new GeometricShapeFactory();
shapeFactory.setNumPoints(parts);
shapeFactory.setCentre(new CoordinateXY(0, 0));
shapeFactory.setSize(10);
return shapeFactory.createCircle();
}
}

Wyświetl plik

@ -0,0 +1,229 @@
package com.onthegomap.planetiler.collection;
import java.util.Arrays;
import java.util.function.IntBinaryOperator;
/**
* A min-heap stored in an array where each element has 4 children.
* <p>
* This is about 5-10% faster than the standard binary min-heap for the case of merging sorted lists.
* <p>
* Ported from <a href=
* "https://github.com/graphhopper/graphhopper/blob/master/core/src/main/java/com/graphhopper/coll/MinHeapWithUpdate.java">GraphHopper</a>
* and:
* <ul>
* <li>modified to use {@code long} values instead of {@code float}</li>
* <li>extracted a common interface for subclass implementations</li>
* <li>modified so that each element has 4 children instead of 2 (improves performance by 5-10%)</li>
* <li>performance improvements to minimize array lookups</li>
* </ul>
*
* @see <a href="https://en.wikipedia.org/wiki/D-ary_heap">d-ary heap (wikipedia)</a>
*/
class ArrayDoubleMinHeap implements DoubleMinHeap {
protected static final int NOT_PRESENT = -1;
protected final int[] posToId;
protected final int[] idToPos;
protected final double[] posToValue;
protected final int max;
protected int size;
private final IntBinaryOperator tieBreaker;
/**
* @param elements the number of elements that can be stored in this heap. Currently the heap cannot be resized or
* shrunk/trimmed after initial creation. elements-1 is the maximum id that can be stored in this heap
*/
ArrayDoubleMinHeap(int elements, IntBinaryOperator tieBreaker) {
// we use an offset of one to make the arithmetic a bit simpler/more efficient, the 0th elements are not used!
posToId = new int[elements + 1];
idToPos = new int[elements + 1];
Arrays.fill(idToPos, NOT_PRESENT);
posToValue = new double[elements + 1];
posToValue[0] = Double.NEGATIVE_INFINITY;
this.max = elements;
this.tieBreaker = tieBreaker;
}
private static int firstChild(int index) {
return (index << 2) - 2;
}
private static int parent(int index) {
return (index + 2) >> 2;
}
@Override
public int size() {
return size;
}
@Override
public boolean isEmpty() {
return size == 0;
}
@Override
public void push(int id, double value) {
checkIdInRange(id);
if (size == max) {
throw new IllegalStateException("Cannot push anymore, the heap is already full. size: " + size);
}
if (contains(id)) {
throw new IllegalStateException("Element with id: " + id +
" was pushed already, you need to use the update method if you want to change its value");
}
size++;
posToId[size] = id;
idToPos[id] = size;
posToValue[size] = value;
percolateUp(size);
}
@Override
public boolean contains(int id) {
checkIdInRange(id);
return idToPos[id] != NOT_PRESENT;
}
@Override
public void update(int id, double value) {
checkIdInRange(id);
int pos = idToPos[id];
if (pos < 0) {
throw new IllegalStateException(
"The heap does not contain: " + id + ". Use the contains method to check this before calling update");
}
double prev = posToValue[pos];
posToValue[pos] = value;
int cmp = compareIdPos(value, prev, id, pos);
if (cmp > 0) {
percolateDown(pos);
} else if (cmp < 0) {
percolateUp(pos);
}
}
@Override
public void updateHead(double value) {
posToValue[1] = value;
percolateDown(1);
}
@Override
public int peekId() {
return posToId[1];
}
@Override
public double peekValue() {
return posToValue[1];
}
@Override
public int poll() {
int id = peekId();
posToId[1] = posToId[size];
posToValue[1] = posToValue[size];
idToPos[posToId[1]] = 1;
idToPos[id] = NOT_PRESENT;
size--;
percolateDown(1);
return id;
}
@Override
public void clear() {
for (int i = 1; i <= size; i++) {
idToPos[posToId[i]] = NOT_PRESENT;
}
size = 0;
}
private void percolateUp(int pos) {
assert pos != 0;
if (pos == 1) {
return;
}
final int id = posToId[pos];
final double val = posToValue[pos];
// the finish condition (index==0) is covered here automatically because we set vals[0]=-inf
int parent;
double parentValue;
while (compareIdPos(val, parentValue = posToValue[parent = parent(pos)], id, parent) < 0) {
posToValue[pos] = parentValue;
idToPos[posToId[pos] = posToId[parent]] = pos;
pos = parent;
}
posToId[pos] = id;
posToValue[pos] = val;
idToPos[posToId[pos]] = pos;
}
private void checkIdInRange(int id) {
if (id < 0 || id >= max) {
throw new IllegalArgumentException("Illegal id: " + id + ", legal range: [0, " + max + "[");
}
}
private void percolateDown(int pos) {
if (size == 0) {
return;
}
assert pos > 0;
assert pos <= size;
final int id = posToId[pos];
final double value = posToValue[pos];
int child;
while ((child = firstChild(pos)) <= size) {
// optimization: this is a very hot code path for performance of k-way merging,
// so manually-unroll the loop over the 4 child elements to find the minimum value
int minChild = child;
double minValue = posToValue[child], childValue;
if (++child <= size) {
if (comparePosPos(childValue = posToValue[child], minValue, child, minChild) < 0) {
minChild = child;
minValue = childValue;
}
if (++child <= size) {
if (comparePosPos(childValue = posToValue[child], minValue, child, minChild) < 0) {
minChild = child;
minValue = childValue;
}
if (++child <= size &&
comparePosPos(childValue = posToValue[child], minValue, child, minChild) < 0) {
minChild = child;
minValue = childValue;
}
}
}
if (compareIdPos(value, minValue, id, minChild) <= 0) {
break;
}
posToValue[pos] = minValue;
idToPos[posToId[pos] = posToId[minChild]] = pos;
pos = minChild;
}
posToId[pos] = id;
posToValue[pos] = value;
idToPos[id] = pos;
}
private int comparePosPos(double val1, double val2, int pos1, int pos2) {
if (val1 < val2) {
return -1;
} else if (val1 == val2 && val1 != Double.NEGATIVE_INFINITY) {
return tieBreaker.applyAsInt(posToId[pos1], posToId[pos2]);
}
return 1;
}
private int compareIdPos(double val1, double val2, int id1, int pos2) {
if (val1 < val2) {
return -1;
} else if (val1 == val2 && val1 != Double.NEGATIVE_INFINITY) {
return tieBreaker.applyAsInt(id1, posToId[pos2]);
}
return 1;
}
}

Wyświetl plik

@ -0,0 +1,81 @@
/*
* Licensed to GraphHopper GmbH under one or more contributor
* license agreements. See the NOTICE file distributed with this work for
* additional information regarding copyright ownership.
*
* GraphHopper GmbH licenses this file to you under the Apache License,
* Version 2.0 (the "License"); you may not use this file except in
* compliance with the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package com.onthegomap.planetiler.collection;
import java.util.function.IntBinaryOperator;
/**
* API for min-heaps that keeps track of {@code int} keys in a range from {@code [0, size)} ordered by {@code long}
* values.
* <p>
* Ported from <a href=
* "https://github.com/graphhopper/graphhopper/blob/master/core/src/main/java/com/graphhopper/coll/MinHeapWithUpdate.java">GraphHopper</a>
* and modified to extract a common interface for subclass implementations.
*/
public interface DoubleMinHeap {
/**
* Returns a new min-heap where each element has 4 children backed by elements in an array.
* <p>
* This is slightly faster than a traditional binary min heap due to a shallower, more cache-friendly memory layout.
*/
static DoubleMinHeap newArrayHeap(int elements, IntBinaryOperator tieBreaker) {
return new ArrayDoubleMinHeap(elements, tieBreaker);
}
int size();
boolean isEmpty();
/**
* Adds an element to the heap, the given id must not exceed the size specified in the constructor. Its illegal to
* push the same id twice (unless it was polled/removed before). To update the value of an id contained in the heap
* use the {@link #update} method.
*/
void push(int id, double value);
/**
* @return true if the heap contains an element with the given id
*/
boolean contains(int id);
/**
* Updates the element with the given id. The complexity of this method is O(log(N)), just like push/poll. Its illegal
* to update elements that are not contained in the heap. Use {@link #contains} to check the existence of an id.
*/
void update(int id, double value);
/**
* Updates the weight of the head element in the heap, pushing it down and bubbling up the new min element if
* necessary.
*/
void updateHead(double value);
/**
* @return the id of the next element to be polled, i.e. the same as calling poll() without removing the element
*/
int peekId();
double peekValue();
/**
* Extracts the element with minimum value from the heap
*/
int poll();
void clear();
}

Wyświetl plik

@ -0,0 +1,302 @@
package com.onthegomap.planetiler.geo;
import com.onthegomap.planetiler.collection.DoubleMinHeap;
import org.locationtech.jts.geom.CoordinateSequence;
import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.geom.LinearRing;
import org.locationtech.jts.geom.Polygon;
import org.locationtech.jts.geom.util.GeometryTransformer;
public class VWSimplifier extends GeometryTransformer {
private double tolerance;
private double k;
public VWSimplifier() {}
public VWSimplifier setTolerance(double tolerance) {
this.tolerance = tolerance;
return this;
}
public VWSimplifier setWeight(double k) {
this.k = k;
return this;
}
private class Vertex {
int heapIdx;
int idx;
double x;
double y;
double area;
Vertex prev;
Vertex next;
Vertex(int idx, CoordinateSequence seq) {
this.idx = idx;
this.x = seq.getX(idx);
this.y = seq.getY(idx);
}
public void remove() {
if (prev != null) {
prev.next = next;
}
if (next != null) {
next.prev = prev;
}
}
public double area() {
if (prev == null || next == null) {
return area = Double.POSITIVE_INFINITY;
}
return area = weight(prev.x, prev.y, x, y, next.x, next.y);
}
}
private static class MinHeap {
Vertex[] array;
int length = 0;
MinHeap(int size) {
array = new Vertex[size + 1];
}
int push(Vertex vertex) {
int idx = length++;
vertex.heapIdx = idx;
array[idx] = vertex;
up(idx);
return length;
}
Vertex pop() {
var removed = array[0];
var object = array[--length];
array[length] = null;
if (length > 0) {
array[object.heapIdx = 0] = object;
down(0);
}
return removed;
}
int remove(Vertex removed) {
int i = removed.heapIdx;
var object = array[--length];
array[length] = null;
if (i != length) {
array[object.heapIdx = i] = object;
if (object.area < removed.area) {
up(i);
} else {
down(i);
}
}
return i;
}
void up(int i) {
var object = array[i];
while (i > 0) {
var up = ((i + 1) >> 1) - 1;
var parent = array[up];
if (object.area >= parent.area) {
break;
}
array[parent.heapIdx = i] = parent;
array[object.heapIdx = i = up] = object;
}
}
void down(int i) {
var object = array[i];
while (true) {
var right = (i + 1) << 1;
var left = right - 1;
var down = i;
var child = array[down];
if (left < length && array[left].area < child.area) {
child = array[down = left];
}
if (right < length && array[right].area < child.area) {
child = array[down = right];
}
if (down == i) {
break;
}
array[child.heapIdx = i] = child;
array[object.heapIdx = i = down] = object;
}
}
}
@Override
protected CoordinateSequence transformCoordinates(CoordinateSequence coords, Geometry parent) {
boolean area = parent instanceof LinearRing;
if (coords.size() == 0) {
return coords;
}
int num = coords.size();
// heap = minHeap(compareArea)
DoubleMinHeap heap = DoubleMinHeap.newArrayHeap(num, Integer::compare);
// maxArea = 0
// double maxArea = 0;
Vertex[] points = new Vertex[num];
// intersecting = []
//
Vertex prev = null;
for (int i = 0; i < num; i++) {
Vertex cur = new Vertex(i, coords);
points[i] = cur;
if (prev != null) {
cur.prev = prev;
prev.next = cur;
heap.push(prev.idx, prev.area());
}
prev = cur;
}
if (prev != null) {
heap.push(prev.idx, prev.area());
}
int left = num;
int min = area ? 4 : 2;
while (!heap.isEmpty()) {
var id = heap.poll();
Vertex point = points[id];
// TODO try minheap
// VW(0) 11M 2.8M 743k 373k 3.5k
// VW(0.1) 11M 3.3M 431k 171k 474
// VW(1) 11M 3.2M 300k 130k 479
// VW(20) 11M 1.5M 263k 122k 479
if (point.area > tolerance || left <= min) {
break;
}
// if (point.area < maxArea) {
// point.area = maxArea;
// } else {
// maxArea = point.area;
// }
// // Check that the new segment doesnt intersect with
// // any existing segments, except for the points
// // immediate neighbours.
// if (intersect(heap, point.previous, point.next))
// intersecting.push(point);
// continue
// // Reattempt to process points that previously would
// // have caused intersections when removed.
// while (i = intersecting.pop()) heap.push(i)
point.remove();
left--;
if (point.prev != null) {
heap.update(point.prev.idx, point.prev.area());
}
if (point.next != null) {
heap.update(point.next.idx, point.next.area());
}
}
//
//
// // make sure we include the first and last points even if they are closer than the simplification threshold
// int num = coords.size();
// DoubleMinHeap heap = DoubleMinHeap.newArrayHeap(num, Integer::compare);
// double maxArea = tolerance;
// // double[] weights = new double[num];
// int[] prev = new int[num];
// int[] next = new int[num];
// // System.err.println("start " + num);
// for (int b = 0; b < num; b++) {
// int a = b - 1;
// int c = b + 1;
// double weight = a < 0 || c >= num ? Double.POSITIVE_INFINITY : weight(coords, a, b, c);
// // weights[b] = weight;
// next[b] = c;
// prev[b] = a;
// heap.push(b, weight);
// // System.err.println(" push " + weight);
// }
// while (!heap.isEmpty()) {
// double weight = heap.peekValue();
// int idx = heap.poll();
// // System.err.println(" poll " + weight);
// if (weight > tolerance) {
// break;
// }
// int a = prev[idx];
// int c = next[idx];
// if (a > 0) {
// // System.err.println(" update before=" + weight(coords, prev[a], a, c));
// heap.update(a, weight(coords, prev[a], a, c));
// }
// if (c < num - 1) {
// // update next (remove c, update c=>weight)
// // System.err.println(" update after=" + weight(coords, a, c, next[c]));
// heap.update(c, weight(coords, a, c, next[c]));
// }
// next[prev[idx]] = c;
// prev[next[idx]] = a;
// }
// for polygons, additionally keep at least 2 intermediate points even if they are below simplification threshold
// to avoid collapse.
MutableCoordinateSequence result = new MutableCoordinateSequence();
for (Vertex point = points[0]; point != null; point = point.next) {
// if (weights[i] > tolerance || (area && i < 3)) {
result.forceAddPoint(point.x, point.y);
// }
}
// System.err.println(result.size());
return result;
}
@Override
protected Geometry transformPolygon(Polygon geom, Geometry parent) {
return geom.isEmpty() ? null : super.transformPolygon(geom, parent);
}
@Override
protected Geometry transformLinearRing(LinearRing geom, Geometry parent) {
boolean removeDegenerateRings = parent instanceof Polygon;
Geometry simpResult = super.transformLinearRing(geom, parent);
if (removeDegenerateRings && !(simpResult instanceof LinearRing)) {
return null;
}
return simpResult;
}
private static double triangleArea(double ax, double ay, double bx, double by, double cx, double cy) {
return Math.abs(((ay - cy) * (bx - cx) + (by - cy) * (cx - ax)) / 2);
}
private static double cos(double ax, double ay, double bx, double by, double cx, double cy) {
double den = Math.hypot(bx - ax, by - ay) * Math.hypot(cx - bx, cy - by),
cos = 0;
if (den > 0) {
cos = Math.clamp((ax - bx) * (cx - bx) + (ay - by) * (cy - by) / den, -1, 1);
}
return cos;
}
private double weight(double cos) {
return (-cos) * k + 1;
}
private double weight(double ax, double ay, double bx, double by, double cx, double cy) {
double area = triangleArea(ax, ay, bx, by, cx, cy);
return k == 0 ? area : (area * weight(cos(ax, ay, bx, by, cx, cy)));
}
private double weight(CoordinateSequence coords, int a, int b, int c) {
return weight(
coords.getX(a), coords.getY(a),
coords.getX(b), coords.getY(b),
coords.getX(c), coords.getY(c)
);
}
}

Wyświetl plik

@ -0,0 +1,107 @@
package com.onthegomap.planetiler.geo;
import static com.onthegomap.planetiler.TestUtils.assertSameNormalizedFeature;
import static com.onthegomap.planetiler.TestUtils.newLineString;
import static com.onthegomap.planetiler.TestUtils.newPolygon;
import static com.onthegomap.planetiler.TestUtils.rectangle;
import org.junit.jupiter.api.Test;
import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.geom.util.AffineTransformation;
public class VisvalingamWhyattSimplifierTest {
final int[] rotations = new int[]{0, 45, 90, 180, 270};
private void testSimplify(Geometry in, Geometry expected, double amount) {
for (int rotation : rotations) {
var rotate = AffineTransformation.rotationInstance(Math.PI * rotation / 180);
assertSameNormalizedFeature(
rotate.transform(expected),
new VWSimplifier().setTolerance(amount).setWeight(0).transform(rotate.transform(in))
);
}
}
@Test
void testSimplify2Points() {
testSimplify(newLineString(
0, 0,
10, 10
), newLineString(
0, 0,
10, 10
), 1);
}
@Test
void testRemoveAPoint() {
testSimplify(newLineString(
0, 0,
5, 0.9,
10, 0
), newLineString(
0, 0,
10, 0
), 5);
}
@Test
void testKeepAPoint() {
testSimplify(newLineString(
0, 0,
5, 1.1,
10, 0
), newLineString(
0, 0,
5, 1.1,
10, 0
), 5);
}
@Test
void testPolygonLeaveAPoint() {
testSimplify(
rectangle(0, 10),
newPolygon(
0, 0,
0, 10,
10, 10,
0, 0
),
200
);
}
@Test
void testLine() {
testSimplify(
newLineString(
0, 0,
1, 0.1,
2, 0,
3, 0.1
),
newLineString(
0, 0,
1, 0.1,
2, 0,
3, 0.1
),
0.09
);
testSimplify(
newLineString(
0, 0,
1, 0.1,
2, 0,
3, 0.1
),
newLineString(
0, 0,
3, 0.1
),
0.11
);
}
}