org.opentripplanner.common.geometry.RecursiveGridIsolineBuilder Maven / Gradle / Ivy
Go to download
Show more of this group Show more artifacts with this name
Show all versions of otp Show documentation
Show all versions of otp Show documentation
The OpenTripPlanner multimodal journey planning system
package org.opentripplanner.common.geometry;
import java.util.ArrayDeque;
import java.util.ArrayList;
import java.util.Collections;
import java.util.Comparator;
import java.util.HashMap;
import java.util.HashSet;
import java.util.List;
import java.util.Map;
import java.util.Queue;
import java.util.Set;
import org.slf4j.Logger;
import org.slf4j.LoggerFactory;
import org.locationtech.jts.algorithm.CGAlgorithms;
import org.locationtech.jts.geom.Coordinate;
import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.geom.GeometryFactory;
import org.locationtech.jts.geom.LinearRing;
import org.locationtech.jts.geom.Polygon;
/**
* Compute isoline based on a zFunc and a set of initial coverage points P0={(x,y)} to seed the
* computation.
*
* This assume we have a z = Fz(x,y) function which can gives a z value for any given point (x,y), z
* can be undefined for some region.
*
* There are many tricks to reduce to the minimum the numbers of Fz samplings, explained in the code
* below.
*
* It will compute an isoline for a given z0 value. The isoline is composed of a list of n polygons,
* CW for normal polygons, CCW for "holes". The isoline computation can be called multiple times on
* the same builder for different z0 values: this will reduce the number of Fz sampling as they are
* cached in the builder. Please note that the initial covering points must touch all isolines you
* want to cover.
*
* @author laurent
*/
public class RecursiveGridIsolineBuilder {
public interface ZFunc {
public long z(Coordinate c);
}
public enum Direction {
UP, DOWN, LEFT, RIGHT;
}
/**
* A fast XY-index for tiles, allowing to use maps.
*/
private static class XYIndex {
private int xIndex;
private int yIndex;
private XYIndex(int xIndex, int yIndex) {
this.xIndex = xIndex;
this.yIndex = yIndex;
}
private final XYIndex translated(int dx, int dy) {
return new XYIndex(this.xIndex + dx, this.yIndex + dy);
}
@Override
public final int hashCode() {
return xIndex >> 16 | yIndex;
}
@Override
public final boolean equals(Object other) {
if (other instanceof XYIndex) {
XYIndex otherTileIndex = (XYIndex) other;
return otherTileIndex.xIndex == this.xIndex && otherTileIndex.yIndex == this.yIndex;
}
return false;
}
@Override
public final String toString() {
return "(" + xIndex + "," + yIndex + ")";
}
}
private static class GridDot {
private XYIndex index;
private long z;
private GridEdge up, down, left, right;
private GridDot(XYIndex index, long z) {
this.index = index;
this.z = z;
}
@Override
public final int hashCode() {
return index.hashCode();
}
@Override
public final boolean equals(Object other) {
if (other instanceof GridDot) {
// Only compare position, not Z value which should be identical
return this.index.equals(((GridDot) other).index);
}
return false;
}
@Override
public final String toString() {
return "[Dot" + index + "," + z + "]";
}
}
private static class GridEdge {
// For horizontal edges, A is left and B is right
// For vertical edges, A is bottom and B is top
private GridDot A, B;
boolean horizontal;
int size;
boolean used = false;
private GridEdge(GridDot A, GridDot B, int size, boolean horizontal) {
// We accept A and B in any order. They must be correctly placed however
if (horizontal && A.index.xIndex > B.index.xIndex || !horizontal
&& A.index.yIndex > B.index.yIndex) {
// Must swap A and B
this.A = B;
this.B = A;
} else {
this.A = A;
this.B = B;
}
this.size = size;
this.horizontal = horizontal;
if (horizontal) {
if (this.A.index.yIndex != this.B.index.yIndex)
throw new AssertionError(
"Building horizontal edge with non horizontally-aligned dots");
if (this.B.index.xIndex - this.A.index.xIndex != size)
throw new AssertionError(
"Building horizontal edge with incorrect size vs dot spacing");
} else {
if (this.A.index.xIndex != this.B.index.xIndex)
throw new AssertionError(
"Building vertical edge with non vertically-aligned dots");
if (this.B.index.yIndex - this.A.index.yIndex != size)
throw new AssertionError(
"Building vertical edge with incorrect size vs dot spacing");
}
}
private final void indexEndPoints() {
if (size != 1)
throw new AssertionError("Can't dot-index edge with size != 1");
if (horizontal) {
A.right = this;
B.left = this;
} else {
A.up = this;
B.down = this;
}
}
private final int cut(long z0) {
if (A.z < z0 && z0 <= B.z)
return 1;
if (B.z < z0 && z0 <= A.z)
return -1;
return 0;
}
@Override
public final int hashCode() {
// Normally size does not matter as we won't
// keep edges from different sizes in the same set/map
// horizontality matter, though.
return A.hashCode() + size + (horizontal ? 0x8000 : 0);
}
@Override
public final boolean equals(Object other) {
if (other instanceof GridEdge) {
// Only compare A position, size and horizontality
GridEdge otherEdge = (GridEdge) other;
return A.equals(otherEdge.A) && size == otherEdge.size
&& horizontal == otherEdge.horizontal;
}
return false;
}
@Override
public final String toString() {
return "[" + (horizontal ? "H" : "V") + "-Edge" + A + "-" + B + " L=" + size + "]";
}
}
private static final Logger LOG = LoggerFactory.getLogger(RecursiveGridIsolineBuilder.class);
/**
* Size (in index-unit) of initial (seed) sampling. This *MUST* be a power of 2. Good values are
* 4 (slower but miss less small islands) or 8 (faster but miss more small islands).
*/
private static final int SIZE_0 = 4;
private ZFunc fz;
private int fzInterpolateCount;
private double dX, dY;
private Coordinate center;
private Map allDots;
private Set initialDots;
private List initialEdges;
private boolean debugSeedGrid = false;
private boolean debugCrossingEdges = false;
private Geometry debugGeometry = null;
// private List p0List;
/**
* Create an object to compute isochrones. One may call several time isochronify on the same
* IsoChronificator object, this will re-use the z = f(x,y) sampling if possible, as they are
* kept in cache.
*
* @param center Center point (eg origin)
* @param fz Function returning the z-value for a xy-coordinate
* @param p0List Initial set of coverage points to seed the heuristics
*/
public RecursiveGridIsolineBuilder( double dX, double dY, Coordinate center, ZFunc fz,
List p0List) {
this.dX = dX;
this.dY = dY;
/*
* Center position only purpose is to serve as a reference value to the XY integer indexing,
* so it only needs to be not too far off to prevent int indexes from overflowing.
*/
this.center = center;
// this.p0List = p0List;
LOG.debug("Center={} dX={} dY={}", this.center, dX, dY);
this.fz = fz;
/* Step 1. SEED (1). Compute initial set of dots. */
allDots = new HashMap(p0List.size() / 2);
initialDots = new HashSet(p0List.size() / 2);
for (Coordinate p0 : p0List) {
XYIndex index0 = getIndex(p0, SIZE_0);
for (int dx = -SIZE_0; dx <= SIZE_0; dx += SIZE_0) {
for (int dy = -SIZE_0; dy <= SIZE_0; dy += SIZE_0) {
// This will always create initial dots
GridDot A = getOrCreateDot(index0.translated(dx, dy));
initialDots.add(A);
}
}
}
/*
* Step 2. SEED (2). Create initial edges from initial dots. There will be slightly less
* edges than dots.
*/
initialEdges = new ArrayList(initialDots.size());
for (GridDot A : initialDots) {
// Horizontal
GridDot B = allDots.get(A.index.translated(SIZE_0, 0));
if (B != null) {
GridEdge e = new GridEdge(A, B, SIZE_0, true);
initialEdges.add(e);
}
// Vertical
B = allDots.get(A.index.translated(0, SIZE_0));
if (B != null) {
GridEdge e = new GridEdge(A, B, SIZE_0, false);
initialEdges.add(e);
}
}
LOG.debug("Created {} dots and {} edges out of {} initial points.", initialDots.size(),
initialEdges.size(), p0List.size());
}
public void setDebugSeedGrid(boolean debugSeedGrid) {
this.debugSeedGrid = debugSeedGrid;
}
public void setDebugCrossingEdges(boolean debugCrossingEdges) {
this.debugCrossingEdges = debugCrossingEdges;
}
public Geometry computeIsoline(long z0) {
fzInterpolateCount = 0;
GeometryFactory geomFactory = new GeometryFactory();
/*
* Step 3. DIVIDE. While there are cutting edges from size > 1 to divide, divide them in
* two. Note: edgesToExpand only contains cutting edges.
*/
Queue edgesToDivide = new ArrayDeque();
edgesToDivide.addAll(initialEdges);
Queue edgesToExpand = new ArrayDeque();
while (!edgesToDivide.isEmpty()) {
GridEdge e = edgesToDivide.remove();
if (e.cut(z0) == 0)
continue;
int size2 = e.size / 2;
XYIndex iC = e.horizontal ? e.A.index.translated(size2, 0) : e.A.index.translated(0,
size2);
GridDot C = getOrCreateDot(iC);
GridEdge e1 = new GridEdge(e.A, C, size2, e.horizontal);
GridEdge e2 = new GridEdge(C, e.B, size2, e.horizontal);
GridEdge eCut = e1.cut(z0) != 0 ? e1 : e2;
if (eCut.cut(z0) == 0)
throw new AssertionError("Edge MUST cut!");
if (size2 == 1) {
edgesToExpand.add(eCut);
} else {
edgesToDivide.add(eCut);
}
}
/*
* Step 4. EXPAND. While there are unprocessed cutting edges, expand them by looking around
* them for cutting edges.
*/
Set finalEdges = new HashSet();
Set finalNonCuttingEdges = new HashSet();
while (!edgesToExpand.isEmpty()) {
GridEdge e = edgesToExpand.remove();
if (finalEdges.add(e) == false) {
// Since we may have duplicates in the toExpand Q,
// we prune-out early processed elements.
continue;
}
/**
* Build the 6 remaining edges (e1..e6) around the 2 squares ABDC and DBFE touching e
*
*
* [Horizontal e] [Vertical e]
*
* C--(e2)--D D--(e3)--B--(e6)--F
* | | | | |
* (e1) (e3) (e2) (e) (e5)
* | | | | |
* A--(e)---B C--(e1)--A--(e4)--E
* | |
* (e4) (e6)
* | |
* E--(e5)--F
*
*/
XYIndex iC = e.horizontal ? e.A.index.translated(0, 1) : e.A.index.translated(-1, 0);
XYIndex iD = e.horizontal ? e.B.index.translated(0, 1) : e.B.index.translated(-1, 0);
XYIndex iE = e.horizontal ? e.A.index.translated(0, -1) : e.A.index.translated(1, 0);
XYIndex iF = e.horizontal ? e.B.index.translated(0, -1) : e.B.index.translated(1, 0);
GridDot C = getOrCreateDot(iC);
GridDot D = getOrCreateDot(iD);
GridDot E = getOrCreateDot(iE);
GridDot F = getOrCreateDot(iF);
GridEdge[] e2List = new GridEdge[] { new GridEdge(e.A, C, 1, !e.horizontal),
new GridEdge(C, D, 1, e.horizontal), new GridEdge(e.B, D, 1, !e.horizontal),
new GridEdge(e.A, E, 1, !e.horizontal), new GridEdge(E, F, 1, e.horizontal),
new GridEdge(e.B, F, 1, !e.horizontal) };
for (GridEdge e2 : e2List) {
if (e2.cut(z0) == 0) {
finalNonCuttingEdges.add(e2);
continue; // Do not cut, OK
}
if (finalEdges.contains(e2))
continue; // Already processed
edgesToExpand.add(e2);
}
}
/*
* Note: Here finalEdges only contains cutting edges, and should end up containing all
* cutting edges that are discoverable.
*/
for (GridEdge e : finalEdges) {
e.indexEndPoints();
}
for (GridEdge e : finalNonCuttingEdges) {
e.indexEndPoints();
}
// For later logs.
int finalEdgesSize = finalEdges.size();
int finalNonCuttingEdgesSize = finalNonCuttingEdges.size();
/*
* Step 5. BUILD. Build polygons from finalEdges set.
*/
List debugGeom = new ArrayList();
if (debugSeedGrid) {
for (GridEdge e : initialEdges) {
Coordinate A = getCoordinate(e.A.index);
Coordinate B = getCoordinate(e.B.index);
debugGeom.add(geomFactory.createLineString(new Coordinate[] { A, B }));
}
// for (Coordinate p0 : p0List) {
// debugGeom.add(geomFactory.createPoint(p0));
// }
}
if (debugCrossingEdges) {
for (GridEdge e : finalEdges) {
Coordinate A = getCoordinate(e.A.index);
Coordinate B = getCoordinate(e.B.index);
Coordinate C = interpolate(A, B, e.A.z, e.B.z, z0);
Coordinate C1 = new Coordinate(C.x + (B.y - A.y) * 0.1, C.y + (B.x - A.x) * 0.1);
Coordinate C2 = new Coordinate(C.x - (B.y - A.y) * 0.1, C.y - (B.x - A.x) * 0.1);
debugGeom
.add(geomFactory.createLineString(new Coordinate[] { A, C, C1, C2, C, B }));
}
}
List retval = new ArrayList();
List rings = new ArrayList();
while (!finalEdges.isEmpty()) {
GridEdge e0 = finalEdges.iterator().next();
List polyPoints = new ArrayList();
int cut = e0.cut(z0);
Direction direction = e0.horizontal ? (cut > 0 ? Direction.UP : Direction.DOWN)
: (cut > 0 ? Direction.LEFT : Direction.RIGHT);
GridEdge e = e0;
while (true) {
// Add a point to polyline
Coordinate cA = getCoordinate(e.A.index);
Coordinate cB = getCoordinate(e.B.index);
Coordinate cC = interpolate(cA, cB, e.A.z, e.B.z, z0);
polyPoints.add(cC);
e.used = true;
finalEdges.remove(e);
// Compute next edge from adjacent tile
// Here e1 is always on e.A side, e2 on e.B side
// and e3 on opposite side of tile from e.
GridEdge e1, e2, e3;
Direction d1, d2; // d3: same direction.
switch (direction) {
default: // Never happen
case UP:
e1 = e.A.up;
d1 = Direction.LEFT;
e2 = e.B.up;
d2 = Direction.RIGHT;
e3 = e1.B.right;
break;
case DOWN:
e1 = e.A.down;
d1 = Direction.LEFT;
e2 = e.B.down;
d2 = Direction.RIGHT;
e3 = e1.A.right;
break;
case LEFT:
e1 = e.A.left;
d1 = Direction.DOWN;
e2 = e.B.left;
d2 = Direction.UP;
e3 = e1.A.up;
break;
case RIGHT:
e1 = e.A.right;
d1 = Direction.DOWN;
e2 = e.B.right;
d2 = Direction.UP;
e3 = e1.B.up;
break;
}
boolean ok1 = e1.cut(z0) != 0 && !e1.used;
boolean ok2 = e2.cut(z0) != 0 && !e2.used;
boolean ok3 = e3.cut(z0) != 0 && !e3.used;
if (ok1 && ok2) {
/*
* We can go either turn, pick the best one from e1 or e2 by looking if C lies
* closer to e.A or e.B. Please note this is approximate only, as we should take
* into account real interpolated position on e1 and e2 to compute segment
* lenght. But this gives good approximated results and is probably sufficient
* given the approximated solution anyway.
*/
double dA = Math.max(Math.abs(cA.x - cC.x), Math.abs(cA.y - cC.y));
double dB = Math.max(Math.abs(cB.x - cC.x), Math.abs(cB.y - cC.y));
if (dA <= dB) {
// C closer to A
e = e1;
direction = d1;
} else {
// C closer to B
e = e2;
direction = d2;
}
} else if (ok1) {
e = e1;
direction = d1;
} else if (ok2) {
e = e2;
direction = d2;
} else if (ok3) {
e = e3;
// Same direction as before
} else {
// This must be the end of the polyline...
break;
}
}
// Close the polyline
polyPoints.add(polyPoints.get(0));
LinearRing ring = geomFactory.createLinearRing(polyPoints
.toArray(new Coordinate[polyPoints.size()]));
rings.add(ring);
}
retval.addAll(punchHoles(geomFactory, rings));
LOG.info("Isochrones: {}+{} Fz samples, {} cutting edges, {} non-cutting edges",
allDots.size(), fzInterpolateCount, finalEdgesSize, finalNonCuttingEdgesSize);
/*
* Step 6. CLEAN. Remove edge index from dots.
*/
for (GridDot A : allDots.values()) {
A.up = A.down = A.right = A.left = null;
}
debugGeometry = geomFactory.createGeometryCollection(debugGeom
.toArray(new Geometry[debugGeom.size()]));
return geomFactory.createGeometryCollection(retval.toArray(new Geometry[retval.size()]));
}
public final Geometry getDebugGeometry() {
return debugGeometry;
}
private final XYIndex getIndex(Coordinate p, int size) {
int xIndex = (int) Math.round((p.x - center.x) / dX);
int yIndex = (int) Math.round((p.y - center.y) / dY);
if (size != 1) {
int size2 = size / 2;
xIndex = (xIndex + (xIndex > 0 ? size2 : -size2)) / size * size;
yIndex = (yIndex + (yIndex > 0 ? size2 : -size2)) / size * size;
}
return new XYIndex(xIndex, yIndex);
}
private final Coordinate getCoordinate(XYIndex index) {
return new Coordinate(index.xIndex * dX + center.x, index.yIndex * dY + center.y);
}
private GridDot getOrCreateDot(XYIndex xyIndex) {
GridDot A = allDots.get(xyIndex);
if (A == null) {
A = new GridDot(xyIndex, fz.z(getCoordinate(xyIndex)));
allDots.put(xyIndex, A);
}
return A;
}
private final Coordinate interpolate(Coordinate A, Coordinate B, long zA, long zB, long z0) {
int n = 0;
while (n < 3 && (zA == Long.MAX_VALUE || zB == Long.MAX_VALUE)) {
Coordinate C = new Coordinate((A.x + B.x) / 2.0, (A.y + B.y) / 2.0);
long zC = fz.z(A);
fzInterpolateCount++;
if (zA == Long.MAX_VALUE && z0 <= zC) {
A = C;
zA = zC;
} else {
B = C;
zB = zC;
}
n++;
}
// Take as fallback position if we are still at +inf at one end z0 * 2
if (zA == Long.MAX_VALUE)
zA = z0 * 2;
if (zB == Long.MAX_VALUE)
zB = z0 * 2;
double k = zB == zA ? 0.5 : (z0 - zA) / (double) (zB - zA);
Coordinate C = new Coordinate(A.x * (1.0 - k) + B.x * k, A.y * (1.0 - k) + B.y * k);
return C;
}
@SuppressWarnings("unchecked")
private final List punchHoles(GeometryFactory geomFactory, List rings) {
List shells = new ArrayList(rings.size());
List holes = new ArrayList(rings.size() / 2);
// 1. Split the polygon list in two: shells and holes (CCW and CW)
for (LinearRing ring : rings) {
if (CGAlgorithms.signedArea(ring.getCoordinateSequence()) > 0.0)
holes.add(ring);
else
shells.add(geomFactory.createPolygon(ring));
}
// 2. Sort the shells based on number of points to optimize step 3.
Collections.sort(shells, new Comparator() {
@Override
public int compare(Polygon o1, Polygon o2) {
return o2.getNumPoints() - o1.getNumPoints();
}
});
for (Polygon shell : shells) {
shell.setUserData(new ArrayList());
}
// 3. For each hole, determine which shell it fits in.
for (LinearRing hole : holes) {
outer: {
// Probably most of the time, the first shell will be the one
for (Polygon shell : shells) {
if (shell.contains(hole)) {
((List) shell.getUserData()).add(hole);
break outer;
}
}
// This should not happen, but do not break bad here
// as loosing a hole is not critical, we still have
// sensible data to return.
LOG.error("Cannot find fitting shell for a hole!");
}
}
// 4. Build the list of punched polygons
List punched = new ArrayList(shells.size());
for (Polygon shell : shells) {
List shellHoles = ((List) shell.getUserData());
punched.add(geomFactory.createPolygon((LinearRing) (shell.getExteriorRing()),
shellHoles.toArray(new LinearRing[shellHoles.size()])));
}
return punched;
}
}