From 6d4f3443261a7fcd74a7bfba0fb506585bfc6569 Mon Sep 17 00:00:00 2001 From: Hannes Janetzek Date: Fri, 24 Jan 2014 07:35:32 +0100 Subject: [PATCH] add douglas-peucker simplifier --- .../org/oscim/utils/geom/GeometryUtils.java | 24 +++ vtm/src/org/oscim/utils/geom/SimplifyDP.java | 144 ++++++++++++++++++ vtm/src/org/oscim/utils/geom/SimplifyVW.java | 99 +++++++++--- 3 files changed, 247 insertions(+), 20 deletions(-) create mode 100644 vtm/src/org/oscim/utils/geom/SimplifyDP.java diff --git a/vtm/src/org/oscim/utils/geom/GeometryUtils.java b/vtm/src/org/oscim/utils/geom/GeometryUtils.java index a7534a80..cdd2c85a 100644 --- a/vtm/src/org/oscim/utils/geom/GeometryUtils.java +++ b/vtm/src/org/oscim/utils/geom/GeometryUtils.java @@ -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; + } } diff --git a/vtm/src/org/oscim/utils/geom/SimplifyDP.java b/vtm/src/org/oscim/utils/geom/SimplifyDP.java new file mode 100644 index 00000000..f8aefee2 --- /dev/null +++ b/vtm/src/org/oscim/utils/geom/SimplifyDP.java @@ -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 . + */ + +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; + } +} diff --git a/vtm/src/org/oscim/utils/geom/SimplifyVW.java b/vtm/src/org/oscim/utils/geom/SimplifyVW.java index bbea6c4e..f6f68e0e 100644 --- a/vtm/src/org/oscim/utils/geom/SimplifyVW.java +++ b/vtm/src/org/oscim/utils/geom/SimplifyVW.java @@ -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 { @@ -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; } } }