Improve LineStringSnapper performance by using squared distance (#1111

)

modules/core/src/main/java/org/locationtech/jts/algorithm/Distance.java

@@ -188,6 +188,64 @@ public static double pointToSegment(Coordinate p, Coordinate A,
         / len2;
     return Math.abs(s) * Math.sqrt(len2);
   }
+
+  /**
+   * Computes the squared distance from a point p to a line segment AB.
+   * 
+   * Note: NON-ROBUST!
+   * 
+   * @param p
+   *          the point to compute the distance for
+   * @param A
+   *          one point of the line
+   * @param B
+   *          another point of the line (must be different to A)
+   * @return the distance from p to line segment AB
+   */
+  public static double pointToSegmentSq(Coordinate p, Coordinate A,
+	      Coordinate B)
+	  {
+	    // if start = end, then just compute distance to one of the endpoints
+	    if (A.x == B.x && A.y == B.y)
+	      return p.distanceSq(A);
+
+	    // otherwise use comp.graphics.algorithms Frequently Asked Questions method
+	    /*
+	     * (1) r = AC dot AB 
+	     *         --------- 
+	     *         ||AB||^2 
+	     *         
+	     * r has the following meaning: 
+	     *   r=0 P = A 
+	     *   r=1 P = B 
+	     *   r<0 P is on the backward extension of AB 
+	     *   r>1 P is on the forward extension of AB 
+	     *   0<r<1 P is interior to AB
+	     */
+
+	    double len2 = (B.x - A.x) * (B.x - A.x) + (B.y - A.y) * (B.y - A.y);
+	    double r = ((p.x - A.x) * (B.x - A.x) + (p.y - A.y) * (B.y - A.y))
+	        / len2;
+
+	    if (r <= 0.0)
+	      return p.distanceSq(A);
+	    if (r >= 1.0)
+	      return p.distanceSq(B);
+
+	    /*
+	     * (2) s = (Ay-Cy)(Bx-Ax)-(Ax-Cx)(By-Ay) 
+	     *         ----------------------------- 
+	     *                    L^2
+	     * 
+	     * Then the distance from C to P = |s|*L.
+	     * 
+	     * This is the same calculation as {@link #distancePointLinePerpendicular}.
+	     * Unrolled here for performance.
+	     */
+	    double s = ((A.y - p.y) * (B.x - A.x) - (A.x - p.x) * (B.y - A.y))
+	        / len2;
+	    return s*s*len2;
+	  }
 
   /**
    * Computes the perpendicular distance from a point p to the (infinite) line

modules/core/src/main/java/org/locationtech/jts/geom/Coordinate.java

@@ -438,6 +438,19 @@ public double distance(Coordinate c) {
     double dy = y - c.y;
     return MathUtil.hypot(dx, dy);
   }
+
+  /**
+   * Computes the 2-dimensional squared Euclidean distance to another location.
+   * The Z-ordinate is ignored.
+   * 
+   * @param c a point
+   * @return the 2-dimensional squared Euclidean distance between the locations
+   */
+  public double distanceSq(Coordinate c) {
+    double dx = x - c.x;
+    double dy = y - c.y;
+    return dx * dx + dy * dy;
+  }
 
   /**
    * Computes the 3-dimensional Euclidean distance to another location.

modules/core/src/main/java/org/locationtech/jts/operation/overlay/snap/LineStringSnapper.java

@@ -12,6 +12,7 @@
 
 package org.locationtech.jts.operation.overlay.snap;
 
+import org.locationtech.jts.algorithm.Distance;
 import org.locationtech.jts.geom.Coordinate;
 import org.locationtech.jts.geom.CoordinateList;
 import org.locationtech.jts.geom.LineSegment;
@@ -30,9 +31,9 @@
 public class LineStringSnapper
 {
   private double snapTolerance = 0.0;
+  private double snapToleranceSq = 0.0;
 
   private Coordinate[] srcPts;
-  private LineSegment seg = new LineSegment(); // for reuse during snapping
   private boolean allowSnappingToSourceVertices = false;
   private boolean isClosed = false;
 
@@ -60,6 +61,7 @@ public LineStringSnapper(Coordinate[] srcPts, double snapTolerance)
     this.srcPts = srcPts;
     isClosed = isClosed(srcPts);
     this.snapTolerance = snapTolerance;
+    this.snapToleranceSq = snapTolerance*snapTolerance;
   }
 
   public void setAllowSnappingToSourceVertices(boolean allowSnappingToSourceVertices)
@@ -191,23 +193,23 @@ private int findSegmentIndexToSnap(Coordinate snapPt, CoordinateList srcCoords)
     double minDist = Double.MAX_VALUE;
     int snapIndex = -1;
     for (int i = 0; i < srcCoords.size() - 1; i++) {
-      seg.p0 = (Coordinate) srcCoords.get(i);
-      seg.p1 = (Coordinate) srcCoords.get(i + 1);
+      final Coordinate p0 = srcCoords.get(i);
+      final Coordinate p1 = srcCoords.get(i+1);
 
       /**
        * Check if the snap pt is equal to one of the segment endpoints.
        * 
        * If the snap pt is already in the src list, don't snap at all.
        */
-      if (seg.p0.equals2D(snapPt) || seg.p1.equals2D(snapPt)) {
+      if (p0.equals2D(snapPt) || p1.equals2D(snapPt)) {
         if (allowSnappingToSourceVertices)
           continue;
         else
           return -1;
       }
 
-      double dist = seg.distance(snapPt);
-      if (dist < snapTolerance && dist < minDist) {
+      double dist = Distance.pointToSegmentSq(snapPt, p0, p1);
+      if (dist < snapToleranceSq && dist < minDist) {
         minDist = dist;
         snapIndex = i;
       }