All Downloads are FREE. Search and download functionalities are using the official Maven repository.

org.opentripplanner.common.geometry.DelaunayIsolineBuilder Maven / Gradle / Ivy

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 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 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 implements IsolineBuilder {

    private static final Logger LOG = LoggerFactory.getLogger(DelaunayIsolineBuilder.class);

    private ZMetric zMetric;

    private DelaunayTriangulation triangulation;

    private boolean debug = false;

    private GeometryFactory geometryFactory = new GeometryFactory();

    private List debugGeom = new ArrayList();

    /**
     * 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 triangulation, ZMetric zMetric) {
        this.triangulation = triangulation;
        this.zMetric = zMetric;
    }

    public void setDebug(boolean debug) {
        this.debug = debug;
    }

    @Override
    public Geometry computeIsoline(TZ z0) {
        Queue> processQ = new ArrayDeque>(
                triangulation.edgesCount());
        for (DelaunayEdge e : triangulation.edges()) {
            e.setProcessed(false);
            processQ.add(e);
        }

        if (debug)
            generateDebugGeometry(z0);

        List rings = new ArrayList();
        while (!processQ.isEmpty()) {
            DelaunayEdge 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 polyPoints = new ArrayList();
            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 E1 = e.getEdge1(ccw);
                DelaunayEdge 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 retval = punchHoles(rings);
        return geometryFactory
                .createGeometryCollection(retval.toArray(new Geometry[retval.size()]));
    }

    private final void generateDebugGeometry(TZ z0) {
        debug = false;
        for (DelaunayEdge 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 punchHoles(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(geometryFactory.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.
        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) 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 punched = new ArrayList(shells.size());
        for (Polygon shell : shells) {
            List shellHoles = ((List) shell.getUserData());
            punched.add(geometryFactory.createPolygon((LinearRing) (shell.getExteriorRing()),
                    shellHoles.toArray(new LinearRing[shellHoles.size()])));
        }
        return punched;
    }
}




© 2015 - 2024 Weber Informatics LLC | Privacy Policy