Add SnappingNoder seeding

modules/app/src/main/java/org/locationtech/jtstest/function/SpatialIndexFunctions.java

@@ -29,13 +29,15 @@
 import org.locationtech.jts.index.chain.MonotoneChain;
 import org.locationtech.jts.index.chain.MonotoneChainBuilder;
 import org.locationtech.jts.index.hprtree.HPRtree;
+import org.locationtech.jts.index.kdtree.KdNode;
 import org.locationtech.jts.index.kdtree.KdTree;
 import org.locationtech.jts.index.quadtree.Quadtree;
 import org.locationtech.jts.index.strtree.AbstractNode;
 import org.locationtech.jts.index.strtree.Boundable;
 import org.locationtech.jts.index.strtree.GeometryItemDistance;
 import org.locationtech.jts.index.strtree.STRtree;
 //import org.locationtech.jts.triangulatepoly.VertexSequencePackedRtree;
+import org.locationtech.jts.math.MathUtil;
 
 
 public class SpatialIndexFunctions
@@ -76,6 +78,117 @@ public static Geometry kdTreeQueryRepeated(Geometry pts, Geometry queryEnv, doub
     return pts.getFactory().createMultiPoint(resultCoords);
   }
 
+  public static Geometry kdTreeGraphSeed(Geometry geom) {
+    return kdTreeGraph(geom, buildKdTreeSeed(geom, 0));
+  }
+
+  public static Geometry kdTreeGraph(Geometry geom) {
+    return kdTreeGraph(geom, buildKdTree(geom, 0));
+  }
+
+  private static Geometry kdTreeGraph(Geometry geom, KdTree index) {
+    KdNode root = index.getRoot();
+    List<Geometry> edges = new ArrayList<Geometry>();
+
+    double x = geom.getEnvelopeInternal().centre().getX();
+    double xInc = geom.getEnvelopeInternal().getWidth() / 2;
+    addGraphEdges(root, true, 0, x, xInc, edges, geom.getFactory());
+    return geom.getFactory().buildGeometry(edges);
+  }
+
+  private static void addGraphEdges(KdNode node, 
+      boolean isXLevel, int depth, double x, double xInc,
+      List<Geometry> edges, GeometryFactory factory) {
+    double xInc2 = xInc / 2;
+    KdNode left = node.getLeft();
+    if (left != null) {
+      double xLeft = x - xInc2;
+      Geometry edgeLeft = factory.createLineString(new Coordinate[] {
+          new Coordinate(x, -depth), new Coordinate(xLeft, -depth-1)
+      });
+      edges.add(edgeLeft);
+      addGraphEdges(left, ! isXLevel, depth+1, xLeft, xInc2, edges, factory);
+    }
+    KdNode right = node.getRight();
+    if (right != null) {
+      double xRight = x + xInc2;
+      Geometry edgeRight = factory.createLineString(new Coordinate[] {
+          new Coordinate(x, -depth), new Coordinate(xRight, -depth-1)
+      });
+      edges.add(edgeRight);
+      addGraphEdges(right, ! isXLevel, depth+1, xRight, xInc2, edges, factory);
+    }
+  }
+
+  public static Geometry kdTreeSplits(Geometry geom) {
+    return kdTreeSplits(geom, buildKdTree(geom, 0));
+  }
+
+  public static Geometry kdTreeSplitsSeed(Geometry geom) {
+    return kdTreeSplits(geom, buildKdTreeSeed(geom, 0));
+  }
+
+  private static Geometry kdTreeSplits(Geometry geom, KdTree index) {
+    Envelope extent = geom.getEnvelopeInternal();
+    KdNode root = index.getRoot();
+    List<Geometry> splits = new ArrayList<Geometry>();
+
+    addSplits(root, true, extent, splits, geom.getFactory());
+    return geom.getFactory().buildGeometry(splits);
+  }
+
+  private static void addSplits(KdNode node, boolean isXLevel, Envelope extent, List<Geometry> splits,
+      GeometryFactory factory) {
+    double splitVal = node.splitValue(isXLevel);
+    Geometry splitLine = createSplitLine(extent, splitVal, isXLevel, factory);
+    splits.add(splitLine);
+
+    KdNode left = node.getLeft();
+    if (left != null) {
+      addSplits(left, ! isXLevel, splitExtent(extent, splitVal, isXLevel, true), splits, factory);
+    }
+    KdNode right = node.getRight();
+    if (right != null) {
+      addSplits(right, ! isXLevel, splitExtent(extent, splitVal, isXLevel, false), splits, factory);
+    }
+  }
+
+  private static Envelope splitExtent(Envelope extent, double splitVal, boolean isXLevel, boolean isLeft) {
+    double xMin = extent.getMinX();
+    double yMin = extent.getMinY();
+    double xMax = extent.getMaxX();
+    double yMax = extent.getMaxY();
+    if (isXLevel) {
+      if (isLeft) {
+        xMax = splitVal;
+      }
+      else {
+        xMin = splitVal;
+      }
+    }
+    else {
+      if (isLeft) {
+        yMax = splitVal;
+      }
+      else {
+        yMin = splitVal;
+      }
+    }
+    return new Envelope(xMin, xMax, yMin, yMax);
+  }
+
+  private static Geometry createSplitLine(Envelope extent, double splitVal, 
+      boolean isXLevel, GeometryFactory factory) {
+
+    double x1 = isXLevel ? splitVal : extent.getMinX();
+    double y1 = isXLevel ? extent.getMinY() : splitVal;
+    double x2 = isXLevel ? splitVal : extent.getMaxX();
+    double y2 = isXLevel ? extent.getMaxY() : splitVal;
+
+    Coordinate[] pts = { new Coordinate(x1, y1), new Coordinate(x2, y2) };
+    return factory.createLineString(pts);
+  }
+
   private static KdTree buildKdTree(Geometry geom, double tolerance) {
     final KdTree index = new KdTree(tolerance);
     Coordinate[] pt = geom.getCoordinates();
@@ -85,6 +198,25 @@ private static KdTree buildKdTree(Geometry geom, double tolerance) {
     return index;
   }
 
+  private static KdTree buildKdTreeSeed(Geometry geom, double tolerance) {
+    final KdTree tree = new KdTree(tolerance);
+    Coordinate[] pt = geom.getCoordinates();
+
+    //-- seed the tree with some randomly selected points
+    int numSeed = pt.length / 100;
+    double rand = 0;
+    for (int i = 0; i < numSeed; i++) {
+      rand = MathUtil.quasirandom(rand);
+      int index = (int) (pt.length * rand);
+      tree.insert(pt[index]);
+    }
+    //-- insert all the points
+    for (int i = 0; i < pt.length; i++) {
+      tree.insert(pt[i]);
+    }
+    return tree;
+  }
+
   public static Geometry strTreeBounds(Geometry geoms)
   {
     STRtree index = new STRtree();

modules/core/src/main/java/org/locationtech/jts/index/kdtree/KdTree.java

@@ -116,6 +116,15 @@ public KdTree(double tolerance) {
     this.tolerance = tolerance;
   }
 
+  /**
+   * Gets the root node of this tree.
+   * 
+   * @return the root node of the tree
+   */
+  public KdNode getRoot() {
+    return root;
+  }
+
   /**
    * Tests whether the index contains any items.
    * 
@@ -453,4 +462,5 @@ private int sizeNode(KdNode currentNode) {
     int sizeR = sizeNode(currentNode.getRight());
     return 1 + sizeL + sizeR;
   }
+
 }

modules/core/src/main/java/org/locationtech/jts/math/MathUtil.java

@@ -147,4 +147,51 @@ public static double min(double v1, double v2, double v3, double v4)
     if (v4 < min) min = v4;
     return min;
   }
+
+  /**
+   * The inverse of the Golden Ratio phi.
+   */
+  public static final double PHI_INV = (Math.sqrt(5) - 1.0) / 2.0;
+
+  /**
+   * Generates a quasi-random sequence of numbers in the range [0,1].
+   * They are produced by an additive recurrence with 1/&phi; as the constant.
+   * This produces a low-discrepancy sequence which is more evenly
+   * distribute than random numbers.
+   * <p>
+   * See <a href='https://en.wikipedia.org/wiki/Low-discrepancy_sequence#Additive_recurrence'>Wikipedia: Low-discrepancy Sequences - Additive Recurrence</a>.
+   * <p>
+   * The sequence is initialized by calling it 
+   * with any positive fractional number; 0 works well for most uses.
+   * 
+   * @param curr the current number in the sequence
+   * @return the next value in the sequence
+   */
+  public static double quasirandom( double curr) {
+    return quasirandom(curr, PHI_INV);
+  }
+
+  /**
+   * Generates a quasi-random sequence of numbers in the range [0,1].
+   * They are produced by an additive recurrence with constant &alpha;.
+   * <pre>
+   *     R(&alpha;) :  t<sub>n</sub> = { t<sub>0</sub> + n&alpha; },  n = 1,2,3,...   
+   * </pre>
+   * When &alpha; is irrational this produces a 
+   * <a href='https://en.wikipedia.org/wiki/Low-discrepancy_sequence#Additive_recurrence'>Low discrepancy sequence</a>
+   *  which is more evenly
+   * distributed than random numbers.
+   * <p>
+   * The sequence is initialized by calling it 
+   * with any positive fractional number. 0 works well for most uses.
+   * 
+   * @param curr the current number in the sequence
+   * @param alpha the sequence additive constant
+   * @return the next value in the sequence
+   */
+  public static double quasirandom( double curr, double alpha) {
+    double next = curr + alpha;
+    if (next < 1) return next;
+    return next - Math.floor(next);
+  }
 }

modules/core/src/main/java/org/locationtech/jts/noding/snap/SnappingNoder.java

@@ -17,6 +17,7 @@
 
 import org.locationtech.jts.geom.Coordinate;
 import org.locationtech.jts.geom.CoordinateList;
+import org.locationtech.jts.math.MathUtil;
 import org.locationtech.jts.noding.MCIndexNoder;
 import org.locationtech.jts.noding.NodedSegmentString;
 import org.locationtech.jts.noding.Noder;
@@ -64,15 +65,18 @@ public SnappingNoder(double snapTolerance) {
   }
 
   /**
+   * Gets the noded result.
+   * 
 	 * @return a Collection of NodedSegmentStrings representing the substrings
-	 * 
 	 */
   public Collection getNodedSubstrings()
   {
     return nodedResult;
   }
 
   /**
+   * Computes the noding of a set of {@link SegmentString}s.
+   * 
    * @param inputSegStrings a Collection of SegmentStrings
    */
   public void computeNodes(Collection inputSegStrings)
@@ -82,13 +86,41 @@ public void computeNodes(Collection inputSegStrings)
   }
 
   private List<NodedSegmentString> snapVertices(Collection<SegmentString> segStrings) {
+    //Stopwatch sw = new Stopwatch(); sw.start();
+    seedSnapIndex(segStrings);
+
     List<NodedSegmentString> nodedStrings = new ArrayList<NodedSegmentString>();
     for (SegmentString ss : segStrings) {
       nodedStrings.add( snapVertices(ss) );
     }
+    //System.out.format("Index depth = %d   Time: %s\n", snapIndex.depth(), sw.getTimeString());
     return nodedStrings;
   }
 
+  /**
+   * Seeds the snap index with a small set of vertices 
+   * chosen quasi-randomly using a low-discrepancy sequence.
+   * Seeding the snap index KdTree induces a more balanced tree. 
+   * This prevents monotonic runs of vertices
+   * unbalancing the tree and causing poor query performance.
+   *  
+   * @param segStrings the segStrings to be noded
+   */
+  private void seedSnapIndex(Collection<SegmentString> segStrings) {
+    final int SEED_SIZE_FACTOR = 100;
+
+    for (SegmentString ss : segStrings) {
+      Coordinate[] pts = ss.getCoordinates();
+      int numPtsToLoad = pts.length / SEED_SIZE_FACTOR;
+      double rand = 0.0;
+      for (int i = 0; i < numPtsToLoad; i++) {
+        rand = MathUtil.quasirandom(rand);
+        int index = (int) (pts.length * rand);
+        snapIndex.snap(pts[index]);
+      }
+    }
+  }
+
   private NodedSegmentString snapVertices(SegmentString ss) {
     Coordinate[] snapCoords = snap(ss.getCoordinates());
     return new NodedSegmentString(snapCoords, ss.getData());

modules/core/src/main/java/org/locationtech/jts/noding/snap/SnappingPointIndex.java

@@ -59,8 +59,22 @@ public Coordinate snap(Coordinate p) {
     return node.getCoordinate();
   }
 
+  /**
+   * Gets the snapping tolerance value for the index.
+   * 
+   * @return the snapping tolerance value
+   */
   public double getTolerance() {
     return snapTolerance;
   }
 
+  /**
+   * Computes the depth of the index tree.
+   * 
+   * @return the depth of the index tree
+   */
+  public int depth() {
+    return snapPointIndex.depth();
+  }
+
 }