add douglas-peucker simplifier

This commit is contained in:
Hannes Janetzek 2014-01-24 07:35:32 +01:00
parent 2241f3f27f
commit 6d4f344326
3 changed files with 247 additions and 20 deletions

View File

@ -108,4 +108,28 @@ public final class GeometryUtils {
return dx * dx + dy * dy;
}
public static double distance(float[] p, int a, int b) {
float dx = p[a] - p[b];
float dy = p[a + 1] - p[b + 1];
return Math.sqrt(dx * dx + dy * dy);
}
public static double dotProduct(float[] p, int a, int b, int c) {
double ab = distance(p, a, b);
double bc = distance(p, b, c);
double d = ab * bc;
double dotp = 0;
if (d > 0) {
dotp = ((p[a] - p[b]) * (p[c] - p[b]) +
(p[a + 1] - p[b + 1]) * (p[c + 1] - p[b + 1]))
/ d;
if (dotp > 1)
dotp = 1;
else if (dotp < 0)
dotp = 0;
}
return dotp;
}
}

View File

@ -0,0 +1,144 @@
package org.oscim.utils.geom;
/*
* Copyright 2012, 2013 Hannes Janetzek
*
* This file is part of the OpenScienceMap project (http://www.opensciencemap.org).
*
* This program is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free Software
* Foundation, either version 3 of the License, or (at your option) any later version.
*
* This program is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
* PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public License along with
* this program. If not, see <http://www.gnu.org/licenses/>.
*/
import static org.oscim.utils.geom.GeometryUtils.squareSegmentDistance;
import org.oscim.core.GeometryBuffer;
/**
* Douglas-Peucker line simplification with stack and some bit twiddling.
*
* based on: https://github.com/mourner/simplify-js,
* https://github.com/ekeneijeoma/simplify-java
*/
public class SimplifyDP {
int[] markers = new int[1];
int[] stack = new int[32];
public void simplify(GeometryBuffer geom, float sqTolerance) {
short[] idx = geom.index;
int inPos = 0;
int outPos = 0;
for (int i = 0, n = idx.length; i < n; i++) {
int len = idx[i];
if (len < 0)
break;
if (len < 6) {
inPos += len;
outPos += len;
continue;
}
int cnt = simplify(geom.points, inPos, len, outPos, sqTolerance);
idx[i] = (short) cnt;
outPos += cnt;
inPos += len;
}
}
public int simplify(float[] points, int inPos, int length, int out, float sqTolerance) {
/* cheap int bitset (use length / 32 ints)
* might should use boolean, or BitSet,
* as this is not really a memory hog :) */
int n = (length >> 5) + 1;
if (markers.length < n) {
markers = new int[n];
} else {
for (int i = 0; i < n; i++)
markers[i] = 0;
}
int first = 0;
int last = inPos + length - 2;
int index = 0;
float maxSqDist;
float sqDist;
int sp = 0;
while (true) {
maxSqDist = 0;
for (int i = first + 2; i < last; i += 2) {
sqDist = squareSegmentDistance(points, i, first, last);
if (sqDist > maxSqDist) {
index = i;
maxSqDist = sqDist;
}
}
if (maxSqDist > sqTolerance) {
/* get marker for 'pos' and shift marker relative
* position into it;
* (pos & 0x1F == pos % 32) */
int pos = index >> 1;
markers[pos >> 5] |= 1 << (pos & 0x1F);
if (sp + 4 == stack.length) {
int tmp[] = new int[stack.length + 64];
System.arraycopy(stack, 0, tmp, 0, stack.length);
stack = tmp;
}
stack[sp++] = first;
stack[sp++] = index;
stack[sp++] = index;
stack[sp++] = last;
}
if (sp == 0)
break;
last = stack[--sp];
first = stack[--sp];
}
points[out++] = points[inPos];
points[out++] = points[inPos + 1];
last = inPos + length - 2;
O: for (int i = 0; i < markers.length; i++) {
int marker = markers[i];
/* point position of this marker */
int mPos = (i << 6);
for (int j = 0; j < 32; j++) {
/* check if marker is set */
if ((marker & (1 << j)) != 0) {
int pos = mPos + (j << 1);
if (pos >= last)
break O;
points[out++] = points[pos];
points[out++] = points[pos + 1];
}
}
}
points[out++] = points[last];
points[out++] = points[last + 1];
return out;
}
}

View File

@ -1,10 +1,15 @@
package org.oscim.utils.geom;
import org.oscim.core.GeometryBuffer;
import org.oscim.utils.GeometryUtils;
import org.oscim.utils.pool.Inlist;
import org.oscim.utils.pool.Pool;
/**
* Visvalingam-Wyatt simplification
*
* ported from:
* https://github.com/mbloch/mapshaper/blob/master/src/mapshaper-visvalingam.js
* */
public class SimplifyVW {
class Item extends Inlist<Item> {
@ -40,16 +45,16 @@ public class SimplifyVW {
Item prev = null;
Item first = null;
Item it;
size = 0;
if (heap.length < geom.pointPos >> 1)
heap = new Item[geom.pointPos >> 1];
first = prev = push(0, Float.MAX_VALUE);
for (int i = 2; i < geom.pointPos - 2; i += 2) {
it = push(i, GeometryUtils.area(geom.points, i - 2, i, i + 2));
it = push(i, area(geom.points, i - 2, i, i + 2));
prev.next = it;
it.prev = prev;
@ -72,33 +77,88 @@ public class SimplifyVW {
it.prev.next = it.next;
it.next.prev = it.prev;
if (it.prev != first)
update(geom, it.prev);
if (it.next != first)
update(geom, it.next);
it = pool.release(it);
}
first.prev.next = null;
first.prev = null;
it = first;
float[] points = new float[geom.pointPos];
System.arraycopy(geom.points, 0, points, 0, geom.pointPos);
geom.clear();
geom.startPolygon();
first.prev.next = null;
it = first;
while (it != null) {
geom.addPoint(geom.points[it.id], geom.points[it.id + 1]);
//System.out.println(first.id + " " + first.area);
float x = points[it.id];
float y = points[it.id + 1];
geom.addPoint(x, y);
it = it.next;
}
//if (merged > 0)
// System.out.println("merged: " + merged);
// if (geom.points[0] == geom.points[geom.pointPos - 2]
// && geom.points[1] == geom.points[geom.pointPos - 1]) {
// System.out.println("same points");
// geom.pointPos -= 2;
// geom.index[geom.indexPos] -= 2;
// }
first = pool.release(first);
}
public static float area(float[] a, int p1, int p2, int p3) {
// float ax = a[p1], ay = a[p1+1];
// float bx = a[p2], by = a[p2+1];
// float cx = a[p3], cy = a[p3+1];
//
// float area = ((ax - cx) * (by - cy) - (bx - cx) * (ay - cy));
//
// area = (area < 0 ? -area : area) * 0.5f;
// float ab = distance2D(ax, ay, bx, by),
// bc = distance2D(bx, by, cx, cy),
// theta, dotp;
//
// if (ab == 0 || bc == 0) {
// theta = 0;
// } else {
// dotp = ((ax - bx) * (cx - bx) + (ay - by) * (cy - by)) / ab * bc;
// if (dotp >= 1) {
// theta = 0;
// } else if (dotp <= -1) {
// theta = Math.PI;
// } else {
// theta = Math.acos(dotp); // consider using other formula at small dp
// }
// }
// return theta;
//
// float angle = innerAngle(ax, ay, bx, by, cx, cy),
// weight = angle < 0.5 ? 0.1 : angle < 1 ? 0.3 : 1;
float area = GeometryUtils.area(a, p1, p2, p3);
//area *= (1 + GeometryUtils.squareSegmentDistance(a, p2, p1, p3) * 10);
float dotp = (float) GeometryUtils.dotProduct(a, p1, p2, p3);
if (dotp > 0)
return area * (1 - dotp);
return area;
}
private void update(GeometryBuffer geom, Item it) {
remove(it);
it.area = GeometryUtils.area(geom.points, it.prev.id, it.id, it.next.id);
it.area = area(geom.points, it.prev.id, it.id, it.next.id);
//System.out.println("< " + it.area);
push(it);
}
@ -154,21 +214,21 @@ public class SimplifyVW {
};
private void up(int i) {
Item obj = heap[i];
Item it = heap[i];
while (i > 0) {
int up = ((i + 1) >> 1) - 1;
Item parent = heap[up];
if (obj.area >= parent.area)
if (it.area >= parent.area)
break;
heap[parent.index = i] = parent;
heap[obj.index = i = up] = obj;
heap[it.index = i = up] = it;
}
}
private void down(int i) {
Item obj = heap[i];
Item it = heap[i];
while (true) {
int right = (i + 1) << 1;
int left = right - 1;
@ -179,15 +239,14 @@ public class SimplifyVW {
if (left < size && heap[left].area < child.area)
child = heap[down = left];
if (right < size && heap[right].area < child.area) {
if (right < size && heap[right].area < child.area)
child = heap[down = right];
}
if (down == i)
break;
heap[child.index = i] = child;
heap[obj.index = down] = obj;
i = down;
heap[it.index = i = down] = it;
}
}
}