/* 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 General Public License for more details. You should have received a copy of the GNU General Public License along with this program. If not, see <http://www.gnu.org/licenses/>. */ package org.opentripplanner.common.geometry; import java.util.ArrayDeque; import java.util.ArrayList; import java.util.Collections; import java.util.Comparator; import java.util.List; import java.util.Queue; import org.slf4j.Logger; import org.slf4j.LoggerFactory; import com.vividsolutions.jts.algorithm.CGAlgorithms; import com.vividsolutions.jts.geom.Coordinate; import com.vividsolutions.jts.geom.Geometry; import com.vividsolutions.jts.geom.GeometryFactory; import com.vividsolutions.jts.geom.LinearRing; import com.vividsolutions.jts.geom.Polygon; /** * Compute isoline based on a Delaunay triangulation of z samplings. * * 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 value: this will reduce the number of Fz sampling as they are * cached in the builder, and reduce the number of time the Delaunay triangulation has to be built. * * The algorithm is rather simple: for each edges of the triangulation check if the edge is * "cutting" (ie crossing the z0 plane). Then start for each unprocessed cutting edge using a walk * algorithm, keeping high z0 always one the same side, to build a set of closed polygons. Then * process each polygons to punch holes: a CW polygon is a hole in a larger CCW polygon, a CCW * polygon is an islan (shell). * * @author laurent */ public class DelaunayIsolineBuilder<TZ> implements IsolineBuilder<TZ> { private static final Logger LOG = LoggerFactory.getLogger(DelaunayIsolineBuilder.class); private ZMetric<TZ> zMetric; private DelaunayTriangulation<TZ> triangulation; private boolean debug = false; private GeometryFactory geometryFactory = new GeometryFactory(); private List<Geometry> debugGeom = new ArrayList<Geometry>(); /** * Create an object to compute isolines. One may call several time computeIsoline on the same * object, with different z0 values. * * @param triangulation The triangulation to process. Must be closed (no edge at the border * should intersect). * @param zMetric The Z metric (intersection detection and interpolation method). */ public DelaunayIsolineBuilder(DelaunayTriangulation<TZ> triangulation, ZMetric<TZ> zMetric) { this.triangulation = triangulation; this.zMetric = zMetric; } public void setDebug(boolean debug) { this.debug = debug; } @Override public Geometry computeIsoline(TZ z0) { Queue<DelaunayEdge<TZ>> processQ = new ArrayDeque<DelaunayEdge<TZ>>( triangulation.edgesCount()); for (DelaunayEdge<TZ> e : triangulation.edges()) { e.setProcessed(false); processQ.add(e); } if (debug) generateDebugGeometry(z0); List<LinearRing> rings = new ArrayList<LinearRing>(); while (!processQ.isEmpty()) { DelaunayEdge<TZ> e = processQ.remove(); if (e.isProcessed()) continue; e.setProcessed(true); int cut = zMetric.cut(e.getA().getZ(), e.getB().getZ(), z0); if (cut == 0) { continue; // While, next edge } List<Coordinate> polyPoints = new ArrayList<Coordinate>(); boolean ccw = cut > 0; while (true) { // Add a point to polyline Coordinate cA = e.getA().getCoordinates(); Coordinate cB = e.getB().getCoordinates(); double k = zMetric.interpolate(e.getA().getZ(), e.getB().getZ(), z0); Coordinate cC = new Coordinate(cA.x * (1.0 - k) + cB.x * k, cA.y * (1.0 - k) + cB.y * k); polyPoints.add(cC); e.setProcessed(true); DelaunayEdge<TZ> E1 = e.getEdge1(ccw); DelaunayEdge<TZ> E2 = e.getEdge2(ccw); int cut1 = E1 == null ? 0 : zMetric.cut(E1.getA().getZ(), E1.getB().getZ(), z0); int cut2 = E2 == null ? 0 : zMetric.cut(E2.getA().getZ(), E2.getB().getZ(), z0); boolean ok1 = cut1 != 0 && !E1.isProcessed(); boolean ok2 = cut2 != 0 && !E2.isProcessed(); if (ok1) { e = E1; ccw = cut1 > 0; } else if (ok2) { e = E2; ccw = cut2 > 0; } else { // This must be the end of the polyline... break; } } // Close the polyline polyPoints.add(polyPoints.get(0)); if (polyPoints.size() > 5) { // If the ring is smaller than 4 points do not add it, // that will remove too small islands or holes. LinearRing ring = geometryFactory.createLinearRing(polyPoints .toArray(new Coordinate[polyPoints.size()])); rings.add(ring); } } List<Polygon> retval = punchHoles(rings); return geometryFactory .createGeometryCollection(retval.toArray(new Geometry[retval.size()])); } private final void generateDebugGeometry(TZ z0) { debug = false; for (DelaunayEdge<TZ> e : triangulation.edges()) { Coordinate cA = e.getA().getCoordinates(); Coordinate cB = e.getB().getCoordinates(); debugGeom.add(geometryFactory.createLineString(new Coordinate[] { cA, cB })); if (zMetric.cut(e.getA().getZ(), e.getB().getZ(), z0) != 0) { double k = zMetric.interpolate(e.getA().getZ(), e.getB().getZ(), z0); Coordinate cC = new Coordinate(cA.x * (1.0 - k) + cB.x * k, cA.y * (1.0 - k) + cB.y * k); debugGeom.add(geometryFactory.createPoint(cC)); } } } public final Geometry getDebugGeometry() { return geometryFactory.createGeometryCollection(debugGeom.toArray(new Geometry[debugGeom .size()])); } @SuppressWarnings("unchecked") private final List<Polygon> punchHoles(List<LinearRing> rings) { List<Polygon> shells = new ArrayList<Polygon>(rings.size()); List<LinearRing> holes = new ArrayList<LinearRing>(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(geometryFactory.createPolygon(ring)); } // 2. Sort the shells based on number of points to optimize step 3. Collections.sort(shells, new Comparator<Polygon>() { @Override public int compare(Polygon o1, Polygon o2) { return o2.getNumPoints() - o1.getNumPoints(); } }); for (Polygon shell : shells) { shell.setUserData(new ArrayList<LinearRing>()); } // 3. For each hole, determine which shell it fits in. int nHolesFailed = 0; 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<LinearRing>) 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. nHolesFailed += 1; } } if (nHolesFailed > 0) { LOG.error("Could not find a shell for {} holes.", nHolesFailed); } // 4. Build the list of punched polygons List<Polygon> punched = new ArrayList<Polygon>(shells.size()); for (Polygon shell : shells) { List<LinearRing> shellHoles = ((List<LinearRing>) shell.getUserData()); punched.add(geometryFactory.createPolygon((LinearRing) (shell.getExteriorRing()), shellHoles.toArray(new LinearRing[shellHoles.size()]))); } return punched; } }