Merge pull request #854 from GIScience/contour_overhaul

Overhaul Contour creation

CHANGELOG.md

@@ -34,6 +34,8 @@ RELEASING:
  -->
 
 ## [Unreleased]
+### Changed
+- Overhaul of Contour creation for fast isochrones. Fixing unexpected behaviour for border edges.
 
 ## [6.3.3] - 2021-01-15
 ### Fixed

openrouteservice/src/main/java/org/heigit/ors/fastisochrones/AbstractEccentricity.java

@@ -1,12 +1,13 @@
 package org.heigit.ors.fastisochrones;
 
+import com.graphhopper.routing.util.EdgeFilter;
 import com.graphhopper.routing.util.FlagEncoder;
 import com.graphhopper.routing.weighting.Weighting;
 import com.graphhopper.storage.GraphHopperStorage;
-import org.heigit.ors.fastisochrones.storage.BorderNodeDistanceStorage;
 import org.heigit.ors.fastisochrones.partitioning.storage.CellStorage;
-import org.heigit.ors.fastisochrones.storage.EccentricityStorage;
 import org.heigit.ors.fastisochrones.partitioning.storage.IsochroneNodeStorage;
+import org.heigit.ors.fastisochrones.storage.BorderNodeDistanceStorage;
+import org.heigit.ors.fastisochrones.storage.EccentricityStorage;
 
 import java.util.ArrayList;
 import java.util.List;
@@ -24,11 +25,11 @@ public abstract class AbstractEccentricity {
     protected List<EccentricityStorage> eccentricityStorages = new ArrayList<>();
     protected List<BorderNodeDistanceStorage> borderNodeDistanceStorages = new ArrayList<>();
 
-    public AbstractEccentricity(GraphHopperStorage ghStorage) {
+    protected AbstractEccentricity(GraphHopperStorage ghStorage) {
         this.ghStorage = ghStorage;
     }
 
-    public abstract void calcEccentricities(Weighting weighting, FlagEncoder flagEncoder);
+    public abstract void calcEccentricities(Weighting weighting, EdgeFilter additionalEdgeFilter, FlagEncoder flagEncoder);
 
     public EccentricityStorage getEccentricityStorage(Weighting weighting) {
         if (eccentricityStorages.isEmpty())

openrouteservice/src/main/java/org/heigit/ors/fastisochrones/Contour.java

@@ -6,6 +6,7 @@
 import com.carrotsearch.hppc.IntSet;
 import com.carrotsearch.hppc.cursors.IntCursor;
 import com.carrotsearch.hppc.cursors.IntObjectCursor;
+import com.graphhopper.routing.util.DefaultEdgeFilter;
 import com.graphhopper.routing.util.EdgeFilter;
 import com.graphhopper.storage.GraphHopperStorage;
 import com.graphhopper.storage.NodeAccess;
@@ -38,13 +39,14 @@
  */
 public class Contour {
     //Length that a polygon edge has to have in m to be split into smaller subedges so that there are no artifacts from later concave hull calculations with it
-    private static final int MIN_EDGE_LENGTH = 400;
+    private static final int MIN_EDGE_LENGTH = 125;
     private static final int MAX_EDGE_LENGTH = Integer.MAX_VALUE;
     //This means that one supercell can contain at most 2^3 = 8 base cells.
-    private static final int SUPER_CELL_HIERARCHY_LEVEL = 3;
+    private static final int SUPER_CELL_HIERARCHY_LEVEL = 2;
     //This means that one supersupercell can contain at most 2^(3 + 2) = 32 base cells.
     private static final int SUPER_SUPER_CELL_HIERARCHY_LEVEL = 2; // level above super cell level
     private static final double CONCAVE_HULL_THRESHOLD = 0.006;
+    private static final double BUFFER_SIZE = 0.0003;
     protected NodeAccess nodeAccess;
     protected GraphHopperStorage ghStorage;
     private IsochroneNodeStorage isochroneNodeStorage;
@@ -97,7 +99,8 @@ public void calculateContour() {
      */
     private void handleBaseCells() {
         for (IntCursor cellId : isochroneNodeStorage.getCellIds()) {
-            LineString ring = createContour(createCoordinates(cellId.value));
+            List<Coordinate> coordinates = createCoordinates(cellId.value);
+            LineString ring = createContour(coordinates, cellStorage.getNodesOfCell(cellId.value).size() < 1000);
             if (ring == null || ring.getNumPoints() < 2) {
                 cellStorage.setCellContourOrder(cellId.value, new ArrayList<>(), new ArrayList<>());
                 continue;
@@ -128,7 +131,7 @@ private IntObjectMap<IntHashSet> handleSuperCells() {
             //Calculate the concave hull for all super cells and super super cells
             for (IntObjectCursor<IntHashSet> superCell : superCellsToBaseCells) {
                 List<Coordinate> superCellCoordinates = createSuperCellCoordinates(superCell.value);
-                LineString ring = createContour(superCellCoordinates);
+                LineString ring = createContour(superCellCoordinates, false);
                 if (ring == null || ring.getNumPoints() < 2) {
                     cellStorage.setCellContourOrder(superCell.key, new ArrayList<>(), new ArrayList<>());
                     continue;
@@ -183,12 +186,12 @@ private List<Coordinate> createSuperCellCoordinates(IntHashSet superCell) {
         return superCellCoordinates;
     }
 
-    private Geometry concHullOfNodes(List<Coordinate> coordinates) {
-        double defaultVisitorThreshold = 0.0035;
+    private Geometry concHullOfNodes(List<Coordinate> coordinates, boolean useHighDetail) {
+        double defaultVisitorThreshold = useHighDetail ? 0.00005 : 0.0025;
         double defaultSearchWidth = 0.0008;
         double defaulPointWidth = 0.005;
 
-        List<Coordinate> points = new ArrayList<>((int)(1 / 20.0 * coordinates.size()));
+        List<Coordinate> points = new ArrayList<>((int) (1 / 20.0 * coordinates.size()));
         PointItemVisitor visitor = new PointItemVisitor(0, 0, defaultVisitorThreshold);
         Quadtree qtree = new Quadtree();
         Envelope searchEnv = new Envelope();
@@ -338,10 +341,8 @@ private void createSuperCell(IntSet cellIds, IntHashSet visitedCells, IntHashSet
         //Cells should be only part of one supercell
         if (visitedCells.contains(currentCell))
             return;
-        if (isPrimary) {
-            if (cellIds.contains(currentCell) && !isValidBaseCell(cellIds, currentCell))
-                return;
-        }
+        if (isPrimary && cellIds.contains(currentCell) && !isValidBaseCell(cellIds, currentCell))
+            return;
 
         if (!cellIds.contains(currentCell)) {
             createSuperCell(cellIds, visitedCells, superCell, maxId, currentCell << 1, isPrimary);
@@ -361,8 +362,9 @@ private List<Coordinate> createCoordinates(int cellId) {
         IntHashSet cellNodes = cellStorage.getNodesOfCell(cellId);
         int initialSize = cellNodes.size();
         List<Coordinate> coordinates = new ArrayList<>(initialSize);
+        EdgeFilter edgeFilter = DefaultEdgeFilter.allEdges(ghStorage.getEncodingManager().fetchEdgeEncoders().get(0));
 
-        EdgeExplorer explorer = ghStorage.getBaseGraph().createEdgeExplorer(EdgeFilter.ALL_EDGES);
+        EdgeExplorer explorer = ghStorage.getBaseGraph().createEdgeExplorer(edgeFilter);
         EdgeIterator iter;
 
         IntHashSet visitedEdges = new IntHashSet();
@@ -371,25 +373,35 @@ private List<Coordinate> createCoordinates(int cellId) {
             towerCoordinates.add(ghStorage.getNodeAccess().getLat(node.value), ghStorage.getNodeAccess().getLon(node.value));
         }
         addLatLon(towerCoordinates, coordinates);
-
         for (IntCursor node : cellNodes) {
             iter = explorer.setBaseNode(node.value);
             while (iter.next()) {
-                if (visitedEdges.contains(iter.getEdge()))
+                if (visitedEdges.contains(iter.getEdge())
+                        || !cellNodes.contains(iter.getAdjNode())
+                        || !edgeFilter.accept(iter))
                     continue;
                 visitedEdges.add(iter.getEdge());
-                addLatLon(iter.fetchWayGeometry(0), coordinates);
+                splitAndAddLatLon(iter.fetchWayGeometry(3), coordinates, MIN_EDGE_LENGTH, MAX_EDGE_LENGTH);
             }
         }
+        //Remove duplicates
+        coordinates = coordinates.stream()
+                .distinct()
+                .collect(Collectors.toList());
         //Need to sort the coordinates, because they will be added to a search tree
         //The order of insertion changes the search tree coordinates and we want consistency between runs
-        Collections.sort(coordinates);
+        try {
+            Collections.sort(coordinates);
+        } catch (Exception e) {
+            //This happens in less than 1% of the runs and I have not figured out why.
+            //It has no real impact as sorting is only for consistency of outlines, not quality
+        }
         return coordinates;
     }
 
-    private LineString createContour(List<Coordinate> coordinates) {
+    private LineString createContour(List<Coordinate> coordinates, boolean useHighDetail) {
         try {
-            Geometry geom = concHullOfNodes(coordinates);
+            Geometry geom = concHullOfNodes(coordinates, useHighDetail);
             Polygon poly = (Polygon) geom;
             poly.normalize();
             return poly.getExteriorRing();
@@ -456,9 +468,44 @@ private void splitEdge(Point point0, Point point1, List<Double> latitudes, List<
         }
     }
 
+    private void splitAndAddLatLon(PointList newCoordinates, List<Coordinate> existingCoordinates, double minlim, double maxlim) {
+        for (int i = 0; i < newCoordinates.size() - 1; i++) {
+            double lat0 = newCoordinates.getLat(i);
+            double lon0 = newCoordinates.getLon(i);
+            double lat1 = newCoordinates.getLat(i + 1);
+            double lon1 = newCoordinates.getLon(i + 1);
+            double dist = distance(lat0, lat1, lon0, lon1);
+            double dx = (lon0 - lon1);
+            double dy = (lat0 - lat1);
+            double normLength = Math.sqrt((dx * dx) + (dy * dy));
+
+            int n = (int) Math.ceil(dist / minlim);
+            double scale = BUFFER_SIZE / normLength;
+
+            double dx2 = -dy * scale;
+            double dy2 = dx * scale;
+            if (i != 0) {
+                existingCoordinates.add(new Coordinate(lon0 + dx2, lat0 + dy2));
+                existingCoordinates.add(new Coordinate(lon0 - dx2, lat0 - dy2));
+            }
+
+            if (dist > minlim && dist < maxlim) {
+                for (int j = 1; j < n; j++) {
+                    existingCoordinates.add(new Coordinate(lon0 + j * (lon1 - lon0) / n + dx2, lat0 + j * (lat1 - lat0) / n + dy2));
+                    existingCoordinates.add(new Coordinate(lon0 + j * (lon1 - lon0) / n - dx2, lat0 + j * (lat1 - lat0) / n - dy2));
+                }
+            }
+        }
+    }
+
     private void addLatLon(PointList newCoordinates, List<Coordinate> existingCoordinates) {
-        for (int i = 0; i < newCoordinates.size(); i++) {
-            existingCoordinates.add(new Coordinate(newCoordinates.getLon(i), newCoordinates.getLat(i)));
+        if (newCoordinates.isEmpty())
+            return;
+        for (int i = 0; i < newCoordinates.size() - 1; i++) {
+            double lat0 = newCoordinates.getLat(i);
+            double lon0 = newCoordinates.getLon(i);
+            existingCoordinates.add(new Coordinate(lon0 + BUFFER_SIZE, lat0 + BUFFER_SIZE));
+            existingCoordinates.add(new Coordinate(lon0 - BUFFER_SIZE, lat0 - BUFFER_SIZE));
         }
     }
 }

openrouteservice/src/main/java/org/heigit/ors/fastisochrones/Eccentricity.java

@@ -12,7 +12,8 @@
 import com.graphhopper.storage.GraphHopperStorage;
 import com.graphhopper.storage.SPTEntry;
 import com.graphhopper.storage.index.LocationIndex;
-import org.heigit.ors.fastisochrones.partitioning.storage.*;
+import org.heigit.ors.fastisochrones.partitioning.storage.CellStorage;
+import org.heigit.ors.fastisochrones.partitioning.storage.IsochroneNodeStorage;
 import org.heigit.ors.fastisochrones.storage.BorderNodeDistanceSet;
 import org.heigit.ors.fastisochrones.storage.BorderNodeDistanceStorage;
 import org.heigit.ors.fastisochrones.storage.EccentricityStorage;
@@ -51,7 +52,7 @@ public Eccentricity(GraphHopperStorage graphHopperStorage, LocationIndex locatio
         this.cellStorage = cellStorage;
     }
 
-    public void calcEccentricities(Weighting weighting, FlagEncoder flagEncoder) {
+    public void calcEccentricities(Weighting weighting, EdgeFilter additionalEdgeFilter, FlagEncoder flagEncoder) {
         if (eccentricityStorages == null) {
             eccentricityStorages = new ArrayList<>();
         }
@@ -83,6 +84,7 @@ public void calcEccentricities(Weighting weighting, FlagEncoder flagEncoder) {
                 FixedCellEdgeFilter fixedCellEdgeFilter = new FixedCellEdgeFilter(isochroneNodeStorage, isochroneNodeStorage.getCellId(node), graph.getNodes());
                 edgeFilterSequence.add(defaultEdgeFilter);
                 edgeFilterSequence.add(fixedCellEdgeFilter);
+                edgeFilterSequence.add(additionalEdgeFilter);
                 RangeDijkstra rangeDijkstra = new RangeDijkstra(graph, weighting);
                 rangeDijkstra.setMaxVisitedNodes(getMaxCellNodesNumber() * eccentricityDijkstraLimitFactor);
                 rangeDijkstra.setEdgeFilter(edgeFilterSequence);
@@ -125,7 +127,7 @@ public void calcEccentricities(Weighting weighting, FlagEncoder flagEncoder) {
         eccentricityStorage.flush();
     }
 
-    public void calcBorderNodeDistances(Weighting weighting, FlagEncoder flagEncoder) {
+    public void calcBorderNodeDistances(Weighting weighting, EdgeFilter additionalEdgeFilter, FlagEncoder flagEncoder) {
         if (borderNodeDistanceStorages == null) {
             borderNodeDistanceStorages = new ArrayList<>();
         }
@@ -140,7 +142,7 @@ public void calcBorderNodeDistances(Weighting weighting, FlagEncoder flagEncoder
         for (IntCursor cellId : isochroneNodeStorage.getCellIds()) {
             final int currentCellId = cellId.value;
             cellCount++;
-            completionService.submit(() -> calculateBorderNodeDistances(borderNodeDistanceStorage, currentCellId, weighting, flagEncoder), String.valueOf(currentCellId));
+            completionService.submit(() -> calculateBorderNodeDistances(borderNodeDistanceStorage, additionalEdgeFilter, currentCellId, weighting, flagEncoder), String.valueOf(currentCellId));
         }
 
         threadPool.shutdown();
@@ -157,14 +159,17 @@ public void calcBorderNodeDistances(Weighting weighting, FlagEncoder flagEncoder
         borderNodeDistanceStorage.flush();
     }
 
-    private void calculateBorderNodeDistances(BorderNodeDistanceStorage borderNodeDistanceStorage, int cellId, Weighting weighting, FlagEncoder flagEncoder) {
+    private void calculateBorderNodeDistances(BorderNodeDistanceStorage borderNodeDistanceStorage, EdgeFilter additionalEdgeFilter, int cellId, Weighting weighting, FlagEncoder flagEncoder) {
         int[] cellBorderNodes = getBorderNodesOfCell(cellId, cellStorage, isochroneNodeStorage).toArray();
+        EdgeFilterSequence edgeFilterSequence = new EdgeFilterSequence();
         EdgeFilter defaultEdgeFilter = DefaultEdgeFilter.outEdges(flagEncoder);
+        edgeFilterSequence.add(defaultEdgeFilter);
+        edgeFilterSequence.add(additionalEdgeFilter);
         Graph graph = ghStorage.getBaseGraph();
 
         for (int borderNode : cellBorderNodes) {
             DijkstraOneToManyAlgorithm algorithm = new DijkstraOneToManyAlgorithm(graph, weighting, TraversalMode.NODE_BASED);
-            algorithm.setEdgeFilter(defaultEdgeFilter);
+            algorithm.setEdgeFilter(edgeFilterSequence);
             algorithm.prepare(new int[]{borderNode}, cellBorderNodes);
             algorithm.setMaxVisitedNodes(getMaxCellNodesNumber() * 20);
             SPTEntry[] targets = algorithm.calcPaths(borderNode, cellBorderNodes);

openrouteservice/src/main/java/org/heigit/ors/fastisochrones/FastIsochroneAlgorithm.java

@@ -21,11 +21,10 @@
 import com.graphhopper.routing.weighting.Weighting;
 import com.graphhopper.storage.Graph;
 import com.graphhopper.storage.SPTEntry;
-import com.vividsolutions.jts.geom.Geometry;
-import org.heigit.ors.fastisochrones.storage.BorderNodeDistanceStorage;
 import org.heigit.ors.fastisochrones.partitioning.storage.CellStorage;
-import org.heigit.ors.fastisochrones.storage.EccentricityStorage;
 import org.heigit.ors.fastisochrones.partitioning.storage.IsochroneNodeStorage;
+import org.heigit.ors.fastisochrones.storage.BorderNodeDistanceStorage;
+import org.heigit.ors.fastisochrones.storage.EccentricityStorage;
 import org.heigit.ors.routing.graphhopper.extensions.edgefilters.EdgeFilterSequence;
 
 import java.util.*;
@@ -39,7 +38,6 @@
 public class FastIsochroneAlgorithm extends AbstractIsochroneAlgorithm {
     private static final String NAME = "FastIsochrone";
     protected IntObjectMap<SPTEntry> startCellMap;
-    protected Set<Integer> activeCells;
     protected Set<Integer> activeBorderNodes;
     protected Set<Integer> inactiveBorderNodes;
     protected Set<Integer> fullyReachableCells;
@@ -67,7 +65,6 @@ protected void initCollections(int size) {
     public void init(int from, double isochroneLimit) {
         this.from = from;
         this.isochroneLimit = isochroneLimit;
-        activeCells = new HashSet<>();
         activeBorderNodes = new HashSet<>();
         inactiveBorderNodes = new HashSet<>();
         fullyReachableCells = new HashSet<>();
@@ -163,10 +160,8 @@ private void findFullyReachableCells(IntObjectMap<SPTEntry> entryMap) {
                     && eccentricityStorage.getFullyReachable(baseNode)) {
                 addFullyReachableCell(baseCell);
                 addInactiveBorderNode(baseNode);
-                removeActiveCell(baseCell);
             } else {
                 if (!getFullyReachableCells().contains(baseCell)) {
-                    addActiveCell(baseCell);
                     addActiveBorderNode(baseNode);
                 }
             }
@@ -175,10 +170,10 @@ private void findFullyReachableCells(IntObjectMap<SPTEntry> entryMap) {
 
     /**
      * Consider all active cells that have a percentage of *approximation* of their nodes visited to be fully reachable.
+     *
      * @param approximation factor of approximation. 1 means all nodes must be found, 0 means no nodes have to be found.
      */
     public void approximateActiveCells(double approximation) {
-        Map<Integer, IntObjectMap<SPTEntry>> newActiveCellMaps = new HashMap<>();
         Iterator<Map.Entry<Integer, IntObjectMap<SPTEntry>>> activeCellIterator = getActiveCellMaps().entrySet().iterator();
         while (activeCellIterator.hasNext()) {
             Map.Entry<Integer, IntObjectMap<SPTEntry>> activeCell = activeCellIterator.next();
@@ -193,14 +188,6 @@ private boolean isWithinLimit(SPTEntry sptEntry, int eccentricity) {
         return sptEntry.getWeightOfVisitedPath() + eccentricity <= isochroneLimit;
     }
 
-    protected void addActiveCell(int cellId) {
-        activeCells.add(cellId);
-    }
-
-    protected void removeActiveCell(int cellId) {
-        activeCells.remove(cellId);
-    }
-
     private void addFullyReachableCell(int cellId) {
         fullyReachableCells.add(cellId);
     }
@@ -217,10 +204,6 @@ public Set<Integer> getFullyReachableCells() {
         return fullyReachableCells;
     }
 
-    private Set<Integer> getActiveCells() {
-        return activeCells;
-    }
-
     public IntObjectMap<SPTEntry> getStartCellMap() {
         return startCellMap;
     }

openrouteservice/src/main/java/org/heigit/ors/fastisochrones/partitioning/FastIsochroneParameters.java

@@ -17,7 +17,7 @@ public final class FastIsochroneParameters {
     //PARTITIONING
     //Based on Implementierung eines Algorithmus zur schnellen Berechnung metrik-affiner Isochronen in einem Straßennetzwerk by Stefan Panig, 2019
     private static int maxCellNodesNumber = 5000;
-    private static int minCellNodesNumber = 4;
+    private static int minCellNodesNumber = 1;
     //Factor based on Aaron Schild & Christian Sommer. On Balanced Seperators in Road Networks, Springer
     //International Publishing Switzerland, 2015 and
     //Implementierung eines Algorithmus zur schnellen Berechnung metrik-affiner Isochronen in einem Straßennetzwerk by Stefan Panig, 2019