From 43b81086c91942a4a1e27ef6c439540e6da6af0c Mon Sep 17 00:00:00 2001 From: Martin Davis Date: Fri, 30 Jan 2026 13:51:50 -0800 Subject: [PATCH 01/43] Add DirectedHausdorffDistance --- .../jtstest/function/DistanceFunctions.java | 25 + .../jtstest/function/SelectionFunctions.java | 11 + .../IndexedPointInPolygonsLocator.java | 2 +- .../distance/DirectedHausdorffDistance.java | 449 ++++++++++++++++++ .../distance/DiscreteHausdorffDistance.java | 1 + .../distance/IndexedFacetDistance.java | 47 +- .../DirectedHausdorffDistanceTest.java | 206 ++++++++ .../distance/HausdorffDistanceStressTest.java | 98 ++++ 8 files changed, 835 insertions(+), 4 deletions(-) create mode 100644 modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java create mode 100644 modules/core/src/test/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistanceTest.java create mode 100644 modules/core/src/test/java/test/jts/perf/operation/distance/HausdorffDistanceStressTest.java diff --git a/modules/app/src/main/java/org/locationtech/jtstest/function/DistanceFunctions.java b/modules/app/src/main/java/org/locationtech/jtstest/function/DistanceFunctions.java index 649db20317..af99d66d36 100644 --- a/modules/app/src/main/java/org/locationtech/jtstest/function/DistanceFunctions.java +++ b/modules/app/src/main/java/org/locationtech/jtstest/function/DistanceFunctions.java @@ -13,6 +13,7 @@ import org.locationtech.jts.algorithm.distance.DiscreteFrechetDistance; import org.locationtech.jts.algorithm.distance.DiscreteHausdorffDistance; +import org.locationtech.jts.algorithm.distance.DirectedHausdorffDistance; import org.locationtech.jts.geom.Coordinate; import org.locationtech.jts.geom.Geometry; import org.locationtech.jts.geom.LineString; @@ -117,4 +118,28 @@ public static Geometry nearestPointsIndexedEachB(Geometry a, Geometry b) { return a.getFactory().createMultiLineString(lines); } + + @Metadata(description="Directed Hausdorff distance from A to B, up to tolerance") + public static double ohd(Geometry a, Geometry b, + @Metadata(title="Distance tolerance") + double distTol) + { + return DirectedHausdorffDistance.distance(a, b, distTol); + } + + @Metadata(description="Directed Hausdorff distance from A to B, up to tolerance") + public static Geometry ohdLine(Geometry a, Geometry b, + @Metadata(title="Distance tolerance") + double distTol) + { + return DirectedHausdorffDistance.distanceLine(a, b, distTol); + } + + @Metadata(description="Hausdorff distance from A to B, up to tolerance") + public static Geometry hdLine(Geometry a, Geometry b, + @Metadata(title="Distance tolerance") + double distTol) + { + return DirectedHausdorffDistance.hausdorffDistanceLine(a, b, distTol); + } } diff --git a/modules/app/src/main/java/org/locationtech/jtstest/function/SelectionFunctions.java b/modules/app/src/main/java/org/locationtech/jtstest/function/SelectionFunctions.java index e3c01b6b19..c24f5d136a 100644 --- a/modules/app/src/main/java/org/locationtech/jtstest/function/SelectionFunctions.java +++ b/modules/app/src/main/java/org/locationtech/jtstest/function/SelectionFunctions.java @@ -16,6 +16,7 @@ import java.util.List; import org.locationtech.jts.algorithm.construct.MaximumInscribedCircle; +import org.locationtech.jts.algorithm.distance.DirectedHausdorffDistance; import org.locationtech.jts.geom.Geometry; import org.locationtech.jts.geom.prep.PreparedGeometry; import org.locationtech.jts.geom.prep.PreparedGeometryFactory; @@ -238,6 +239,16 @@ public boolean isTrue(Geometry g) { }); } + public static Geometry fullyWithinDistance(Geometry a, final Geometry mask, double maximumDistance) + { + return select(a, new GeometryPredicate() { + public boolean isTrue(Geometry g) { + double tolerance = maximumDistance / 1e5; + return DirectedHausdorffDistance.isFullyWithinDistance(a, mask, maximumDistance, tolerance); + } + }); + } + public static Geometry maxInCircleRadiusWithin(Geometry a, @Metadata(title="Max Radius Length") double maximumRadius) diff --git a/modules/core/src/main/java/org/locationtech/jts/algorithm/construct/IndexedPointInPolygonsLocator.java b/modules/core/src/main/java/org/locationtech/jts/algorithm/construct/IndexedPointInPolygonsLocator.java index ce99843041..a5745c5102 100644 --- a/modules/core/src/main/java/org/locationtech/jts/algorithm/construct/IndexedPointInPolygonsLocator.java +++ b/modules/core/src/main/java/org/locationtech/jts/algorithm/construct/IndexedPointInPolygonsLocator.java @@ -29,7 +29,7 @@ * @author mdavis * */ -class IndexedPointInPolygonsLocator implements PointOnGeometryLocator { +public class IndexedPointInPolygonsLocator implements PointOnGeometryLocator { private Geometry geom; private STRtree index; diff --git a/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java b/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java new file mode 100644 index 0000000000..a6d25a079b --- /dev/null +++ b/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java @@ -0,0 +1,449 @@ +/* + * Copyright (c) 2026 Martin Davis + * + * All rights reserved. This program and the accompanying materials + * are made available under the terms of the Eclipse Public License 2.0 + * and Eclipse Distribution License v. 1.0 which accompanies this distribution. + * The Eclipse Public License is available at http://www.eclipse.org/legal/epl-v20.html + * and the Eclipse Distribution License is available at + * + * http://www.eclipse.org/org/documents/edl-v10.php. + */ +package org.locationtech.jts.algorithm.distance; + +import java.util.Iterator; +import java.util.PriorityQueue; + +import org.locationtech.jts.algorithm.construct.IndexedPointInPolygonsLocator; +import org.locationtech.jts.algorithm.construct.LargestEmptyCircle; +import org.locationtech.jts.geom.Coordinate; +import org.locationtech.jts.geom.Dimension; +import org.locationtech.jts.geom.Geometry; +import org.locationtech.jts.geom.GeometryCollectionIterator; +import org.locationtech.jts.geom.GeometryComponentFilter; +import org.locationtech.jts.geom.LineString; +import org.locationtech.jts.geom.Location; +import org.locationtech.jts.geom.Point; +import org.locationtech.jts.operation.distance.IndexedFacetDistance; + +/** + * Computes the directed Hausdorff distance from one geometry to another, + * up to an approximation distance tolerance. + * The directed Hausdorff distance (DHD) is defined as: + * 
+ * DHD(A,B) = maxa ∈ A (maxb ∈ B (distance(a, b) )
+ *
+ * The DHD is the maximum distance any point + * on a query geometry A can be from a target geometry B. + * Equivalently, every point in the query geometry is within the DHD distance + * of the target geometry. + * The class can compute a pair of points at which the DHD is obtained: + * [ farthest A point, nearest B point ]. + *

+ * The DHD is asymmetric: DHD(A,B) may not be equal to DHD(B,A). + * Hence it is not a distance metric. + * The Hausdorff distance is is a symmetric distance metric: + * 

+ * HD(A,B) = max(DHD(A,B), DHD(B,A))
+ *
+ * This can be computed via the + * {@link #hausdorffDistanceLine(Geometry, Geometry, double)} function. + *

+ * Points, lines and polygons are supported as input. + * If the query geometry is polygonal, + * the point of maximum distance may occur in the interior of a polygon. + *

+ * The class can be used in prepared mode. + * Creating an instance on a target geometry caches indexes for that geometry. + * Then {@link #computeDistancePoints(Geometry, double) + * or {@link #isFullyWithinDistance(Geometry, double, double)} + * can be called efficiently for multiple query geometries. + *

+ * A use case is to test whether a geometry A lies fully within a given + * distance of another one B. + * Using {@link #isFullyWithinDistance(Geometry, double, double)} + * is much more efficient than computing whether A is covered by a buffer of B. + *

+ * Due to the nature of the Hausdorff distance, + * performance is not very sensitive to the distance tolerance, + * so using a small tolerance is recommended. + *

+ * This algorithm is easier to use, more accurate, + * and much faster than {@link DiscreteHausdorrffDistance}. + * + * @author Martin Davis + * + */ +public class DirectedHausdorffDistance { + + public static double distance(Geometry a, Geometry b, double tolerance) + { + DirectedHausdorffDistance hd = new DirectedHausdorffDistance(b); + return distance(hd.computeDistancePoints(a, tolerance)); + } + + public static LineString distanceLine(Geometry a, Geometry b, double tolerance) + { + DirectedHausdorffDistance hd = new DirectedHausdorffDistance(b); + return a.getFactory().createLineString(hd.computeDistancePoints(a, tolerance)); + } + + public static boolean isFullyWithinDistance(Geometry a, Geometry b, double maxDistance, double tolerance) + { + DirectedHausdorffDistance hd = new DirectedHausdorffDistance(b); + return hd.isFullyWithinDistance(a, maxDistance, tolerance); + } + + public static LineString hausdorffDistanceLine(Geometry a, Geometry b, double tolerance) + { + DirectedHausdorffDistance hdAB = new DirectedHausdorffDistance(b); + Coordinate[] ptsAB = hdAB.computeDistancePoints(a, tolerance); + DirectedHausdorffDistance hdBA = new DirectedHausdorffDistance(a); + Coordinate[] ptsBA = hdBA.computeDistancePoints(b, tolerance); + + //-- return points in A-B order + Coordinate[] pts; + if (distance(ptsAB) > distance(ptsBA)) { + pts = ptsAB; + } + else { + pts = copyReverse(ptsBA); + } + return a.getFactory().createLineString(pts); + } + + private Geometry geomB; + private IndexedFacetDistance distToB; + private boolean isAreaB; + private IndexedPointInPolygonsLocator ptInAreaB; + + public DirectedHausdorffDistance(Geometry b) { + geomB = b; + distToB = new IndexedFacetDistance(b); + isAreaB = b.getDimension() >= Dimension.A; + if (isAreaB) { + ptInAreaB = new IndexedPointInPolygonsLocator(b); + } + } + + public boolean isFullyWithinDistance(Geometry a, double maxDistance, double tolerance) { + //TODO: optimize with short-circuiting + Coordinate[] maxDistCoords = computeDistancePoints(a, tolerance); + return distance(maxDistCoords) <= maxDistance; + } + + private static double distance(Coordinate[] pts) { + return pts[0].distance(pts[1]); + } + + private Coordinate[] computeDistancePoints(Geometry geomA, double tolerance) { + if (geomA.getDimension() == Dimension.P) { + return computeAtPoints(geomA); + } + //TODO: handle mixed geoms with points + Coordinate[] maxDistPtsEdge = computeAtEdges(geomA, tolerance); + + /** + * Polygonal query geometry may have an interior point as the farthest point. + */ + if (geomA.getDimension() == Dimension.A) { + Coordinate[] maxDistPtsInterior = computeAtAreaInterior(geomA, tolerance); + if (maxDistPtsInterior != null + && distance(maxDistPtsInterior) > distance(maxDistPtsEdge)) { + return maxDistPtsInterior; + } + } + return maxDistPtsEdge; + } + + /** + * If the query geometry A is polygonal, it is possible + * the farthest point lies in its interior. + * In this case it occurs at the centre of the Largest Empty Circle + * with B as obstacles and query geometry as constraint. + * + * This is a potentially expensive computation, + * especially for small tolerances. + * It is avoided where possible if heuristic checks + * can determine that the max distance is at + * the previously computed edge points. + * + * @param geomA + * @param tolerance + * @return the maximum distance point pair at an interior point of A, + * or null if it is known to not occur at an interior point + */ + private Coordinate[] computeAtAreaInterior(Geometry geomA, double tolerance) { + //TODO: extract polygonal geoms from A + Geometry polygonalA = geomA; + + /** + * Optimization - skip if A interior cannot intersect B, + * and thus farther point must lie on A segment + */ + if (polygonalA.getEnvelopeInternal().disjoint(geomB.getEnvelopeInternal())) { + return null; + } + + Point centerPt = LargestEmptyCircle.getCenter(geomB, polygonalA, tolerance); + Coordinate ptA = centerPt.getCoordinate(); + /** + * If LEC centre is in B, the max distance is zero, so return null. + * This will cause the computed segment distance to be returned, + * which is better since it occurs on a vertex (or edge?) + */ + if (isInteriorB(ptA)) { + return null; + } + Coordinate[] nearPtsBA = distToB.nearestPoints(ptA); + return copyReverse(nearPtsBA); + } + + private Coordinate[] computeAtPoints(Geometry geomA) { + double maxDist = -1.0;; + Coordinate[] maxDistPtsAB = null; + Iterator geomi = new GeometryCollectionIterator(geomA); + while (geomi.hasNext()) { + Geometry geomElemA = (Geometry) geomi.next(); + if (! (geomElemA instanceof Point)) + continue; + + Coordinate pA = geomElemA.getCoordinate(); + Coordinate[] nearPtsBA = distToB.nearestPoints(pA); + double dist = distance(nearPtsBA); + + boolean isInterior = dist > 0 && isInteriorB(pA); + //-- check for interior point + if (isInterior) { + dist = 0; + nearPtsBA = copyPair( pA, pA ); + } + if (dist > maxDist) { + maxDist = dist; + maxDistPtsAB = copyReverse(nearPtsBA); + } + } + return maxDistPtsAB; + } + + private static Coordinate[] copyPair(Coordinate p0, Coordinate p1) { + return new Coordinate[] { p0.copy(), p1.copy() }; + } + + private static Coordinate[] copyReverse(Coordinate[] pts) { + return new Coordinate[] { pts[1].copy(), pts[0].copy() }; + } + + public Coordinate[] computeAtEdges(Geometry geomA, double tolerance) { + PriorityQueue segQueue = createSegQueue(geomA); + + DHDSegment ohdSeg = null; + DHDSegment maxDistSeg = null; + long iter = 0; + while (! segQueue.isEmpty()) { + iter++; + // get the segment with greatest distance + maxDistSeg = segQueue.remove(); + + if (maxDistSeg.length() <= tolerance) { + ohdSeg = maxDistSeg; + break; + } + + //-- not within tolerance, so bisect segment and keep searching + DHDSegment[] bisects = maxDistSeg.bisect(distToB); + addNonInterior(bisects[0], segQueue); + addNonInterior(bisects[1], segQueue); + } + if (ohdSeg != null) + return ohdSeg.getMaxDistPts(); + + /** + * If no OHD segment was found, all were inside the target. + * In this case distance is zero. + * Return a single coordinate of the input as a representative point + */ + Coordinate maxPt = maxDistSeg.getEndpoint(0); + return copyPair(maxPt, maxPt); + } + + private void addNonInterior(DHDSegment ohdSegment, PriorityQueue segQueue) { + //-- discard segment if it is interior to a polygon + if (isInterior(ohdSegment)) { + return; + } + segQueue.add(ohdSegment); + } + + /** + * Tests if segment is fully in the interior of the target geometry polygons (if any). + * + * @param ohdSegment + * @return + */ + private boolean isInterior(DHDSegment ohdSegment) { + if (! isAreaB) + return false; + double segDist = distToB.distance(ohdSegment.getEndpoint(0), ohdSegment.getEndpoint(1)); + //-- if segment touches B it is not in interior + if (segDist == 0) + return false; + //-- only need to test one point to check interior + Coordinate pt = ohdSegment.getEndpoint(0); + boolean isInterior = isInteriorB(pt); + return isInterior; + } + + private boolean isInteriorB(Coordinate p) { + if (! isAreaB) return false; + return ptInAreaB.locate(p) == Location.INTERIOR; + } + + private PriorityQueue createSegQueue(Geometry geomA) { + PriorityQueue priq = new PriorityQueue(); + geomA.apply(new GeometryComponentFilter() { + + @Override + public void filter(Geometry geom) { + if (geom instanceof LineString) { + initSegments(geom.getCoordinates(), priq); + } + } + }); + return priq; + } + + /** + * + * @param pts + * @param priq + */ + private void initSegments(Coordinate[] pts, PriorityQueue priq) { + DHDSegment prevSeg = null; + for (int i = 0; i < pts.length - 1; i++) { + DHDSegment seg; + if (i == 0) { + seg = DHDSegment.create(pts[i], pts[i + 1], distToB); + } + else { + //-- optimize by avoiding recomputing pt distance + seg = DHDSegment.create(prevSeg, pts[i + 1], distToB); + } + prevSeg = seg; + /** + * Delay interior check until splitting, + * to avoid computing more than needed. + */ + priq.add(seg); + //System.out.println(seg.distance()); + } + } + + private static class DHDSegment implements Comparable { + + public static DHDSegment create(Coordinate p0, Coordinate p1, IndexedFacetDistance indexDist) { + DHDSegment seg = new DHDSegment(p0, p1); + seg.init(indexDist); + return seg; + } + + public static DHDSegment create(DHDSegment prevSeg, Coordinate p1, IndexedFacetDistance indexDist) { + DHDSegment seg = new DHDSegment(prevSeg.p1, p1); + seg.init(prevSeg.nearPt1, indexDist); + return seg; + } + + private Coordinate p0; + private Coordinate nearPt0; + private Coordinate p1; + private Coordinate nearPt1; + double maxDistanceBound = Double.NEGATIVE_INFINITY; + + private DHDSegment(Coordinate p0, Coordinate p1) { + this.p0 = p0; + this.p1 = p1; + } + + private DHDSegment(Coordinate p0, Coordinate nearPt0, Coordinate p1, Coordinate nearPt1) { + this.p0 = p0; + this.nearPt0 = nearPt0; + this.p1 = p1; + this.nearPt1 = nearPt1; + computeMaxDistanceBound(); + } + + private void init(IndexedFacetDistance indexDist) { + nearPt0 = nearestPt(p0, indexDist); + nearPt1 = nearestPt(p1, indexDist); + computeMaxDistanceBound(); + } + + private void init(Coordinate nearest0, IndexedFacetDistance indexDist) { + nearPt0 = nearest0; + nearPt1 = nearestPt(p1, indexDist); + computeMaxDistanceBound(); + } + + private static Coordinate nearestPt(Coordinate p, IndexedFacetDistance indexDist) { + Coordinate[] nearestPts = indexDist.nearestPoints(p); + return nearestPts[0]; + } + + public Coordinate getEndpoint(int index) { + return index == 0 ? p0 : p1; + } + + public Coordinate getNearestPt(int index) { + return index == 0 ? nearPt0 : nearPt1; + } + + public double length() { + return p0.distance(p1); + } + + public Coordinate[] getMaxDistPts() { + double dist0 = p0.distance(nearPt0); + double dist1 = p1.distance(nearPt1); + if (dist0 > dist1) { + return copyPair(p0, nearPt0); + } + return copyPair(p1, nearPt1); + } + + /** + * Computes a least upper bound for the maximum distance to a segment. + */ + private void computeMaxDistanceBound() { + //System.out.println(distance()); + + /** + * Least upper bound is the max distance to the endpoints, + * plus half segment length. + */ + double dist0 = p0.distance(nearPt0); + double dist1 = p1.distance(nearPt1); + double maxDist = Math.max(dist0, dist1); + maxDistanceBound = maxDist + length() / 2; + } + + public DHDSegment[] bisect(IndexedFacetDistance distToB) { + Coordinate mid = new Coordinate( + (p0.x + p1.x) / 2, + (p0.y + p1.y) / 2 + ); + Coordinate nearPtMid = nearestPt(mid, distToB); + return new DHDSegment[] { + new DHDSegment(p0, nearPt0, mid, nearPtMid ), + new DHDSegment(mid, nearPtMid, p1, nearPt1) + }; + } + + /** + * For maximum efficiency sort the PriorityQueue with largest maxDistance at front. + * Since Java PQ sorts least-first, need to invert the comparison + */ + public int compareTo(DHDSegment o) { + return -Double.compare(maxDistanceBound, o.maxDistanceBound); + } + } +} diff --git a/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DiscreteHausdorffDistance.java b/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DiscreteHausdorffDistance.java index 46a71801cf..b39d923e1f 100644 --- a/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DiscreteHausdorffDistance.java +++ b/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DiscreteHausdorffDistance.java @@ -71,6 +71,7 @@ * This is more efficient than testing whether A is covered by a buffer of B. * * @see DiscreteFrechetDistance + * @see DirectedHausdorffDistance * */ public class DiscreteHausdorffDistance diff --git a/modules/core/src/main/java/org/locationtech/jts/operation/distance/IndexedFacetDistance.java b/modules/core/src/main/java/org/locationtech/jts/operation/distance/IndexedFacetDistance.java index e212ef4785..4ee62e5dc9 100644 --- a/modules/core/src/main/java/org/locationtech/jts/operation/distance/IndexedFacetDistance.java +++ b/modules/core/src/main/java/org/locationtech/jts/operation/distance/IndexedFacetDistance.java @@ -13,10 +13,12 @@ package org.locationtech.jts.operation.distance; import org.locationtech.jts.geom.Coordinate; +import org.locationtech.jts.geom.CoordinateSequence; import org.locationtech.jts.geom.Geometry; import org.locationtech.jts.geom.Lineal; import org.locationtech.jts.geom.Polygonal; import org.locationtech.jts.geom.Puntal; +import org.locationtech.jts.geom.impl.CoordinateArraySequence; import org.locationtech.jts.index.strtree.ItemBoundable; import org.locationtech.jts.index.strtree.ItemDistance; import org.locationtech.jts.index.strtree.STRtree; @@ -87,7 +89,7 @@ public static boolean isWithinDistance(Geometry g1, Geometry g2, double distance * * @param g1 a geometry * @param g2 a geometry - * @return the nearest points on the facets of the geometries + * @return the nearest points on the facets of geometry g1 and g2 */ public static Coordinate[] nearestPoints(Geometry g1, Geometry g2) { IndexedFacetDistance dist = new IndexedFacetDistance(g1); @@ -140,10 +142,10 @@ private Object[] nearestFacets(Geometry g) { /** * Computes the nearest points on the base geometry - * and the given geometry. + * and another geometry. * * @param g the geometry to compute the nearest location to - * @return the nearest points + * @return the nearest points on the base and the argument geometry */ public Coordinate[] nearestPoints(Geometry g) { @@ -153,6 +155,44 @@ public Coordinate[] nearestPoints(Geometry g) return fs1.nearestLocations(fs2); } + /** + * Computes the nearest points on the base geometry + * and a point. + * + * @param p the point coordinate + * @return the nearest point on the base geometry and the point + */ + public Coordinate[] nearestPoints(Coordinate p) + { + CoordinateSequence seq = new CoordinateArraySequence(new Coordinate[] { p }); + FacetSequence fs2 = new FacetSequence(seq, 0); + Object nearest = cachedTree.nearestNeighbour(fs2.getEnvelope(), fs2, FACET_SEQ_DIST); + FacetSequence fs1 = (FacetSequence) nearest; + return fs1.nearestLocations(fs2); + } + + + public double distance(Coordinate p) { + return distance(nearestPoints(p)); + } + + private static double distance(Coordinate[] pts) { + return pts[0].distance(pts[1]); + } + + public Coordinate[] nearestPoints(Coordinate p0, Coordinate p1) + { + CoordinateSequence seq = new CoordinateArraySequence(new Coordinate[] { p0, p1 }); + FacetSequence fs2 = new FacetSequence(seq, 0, 2); + Object nearest = cachedTree.nearestNeighbour(fs2.getEnvelope(), fs2, FACET_SEQ_DIST); + FacetSequence fs1 = (FacetSequence) nearest; + return fs1.nearestLocations(fs2); + } + + public double distance(Coordinate p0, Coordinate p1) { + return distance(nearestPoints(p0, p1)); + } + /** * Tests whether the base geometry lies within * a specified distance of the given geometry. @@ -181,6 +221,7 @@ public double distance(ItemBoundable item1, ItemBoundable item2) { return fs1.distance(fs2); } } + } diff --git a/modules/core/src/test/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistanceTest.java b/modules/core/src/test/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistanceTest.java new file mode 100644 index 0000000000..5d17932b33 --- /dev/null +++ b/modules/core/src/test/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistanceTest.java @@ -0,0 +1,206 @@ +/* + * Copyright (c) 2026 Martin Davis. + * + * All rights reserved. This program and the accompanying materials + * are made available under the terms of the Eclipse Public License 2.0 + * and Eclipse Distribution License v. 1.0 which accompanies this distribution. + * The Eclipse Public License is available at http://www.eclipse.org/legal/epl-v20.html + * and the Eclipse Distribution License is available at + * + * http://www.eclipse.org/org/documents/edl-v10.php. + */ +package org.locationtech.jts.algorithm.distance; + +import org.locationtech.jts.geom.Coordinate; +import org.locationtech.jts.geom.Geometry; + +import junit.textui.TestRunner; +import test.jts.GeometryTestCase; + +public class DirectedHausdorffDistanceTest +extends GeometryTestCase +{ + public static void main(String args[]) { + TestRunner.run(DirectedHausdorffDistanceTest.class); + } + + public DirectedHausdorffDistanceTest(String name) { super(name); } + + public void testPoints() + { + checkHD("POINT (0 0)", "POINT (1 1)", + 0.01, + "LINESTRING (0 0, 1 1)"); + } + + public void testMultiPoints() + { + String a = "MULTIPOINT ((0 1), (2 3), (4 5), (6 6))"; + String b = "MULTIPOINT ((0.1 0), (1 0), (2 0), (3 0), (4 0), (5 0))"; + checkDistance(a, b, 0.01, "LINESTRING (6 6, 5 0)"); + checkDistance(b, a, 0.01, "LINESTRING (5 0, 2 3)"); + checkHD(a, b, 0.01, "LINESTRING (6 6, 5 0)"); + } + + public void testLineSegments() + { + checkHD("LINESTRING (0 0, 2 0)", "LINESTRING (0 0, 2 1)", + 0.01, + "LINESTRING (2 0, 2 1)"); + } + + public void testLineSegments2() + { + checkHD("LINESTRING (0 0, 2 0)", "LINESTRING (0 1, 1 2, 2 1)", + 0.01, + "LINESTRING (1 0, 1 2)"); + } + + public void testLinePoints() + { + checkHD("LINESTRING (0 0, 2 0)", "MULTIPOINT (0 2, 1 0, 2 1)", + 0.01, + "LINESTRING (0 0, 0 2)"); + } + + public void testLinesPolygon() + { + checkHD("MULTILINESTRING ((1 1, 2 7), (7 1, 9 9))", + "POLYGON ((3 7, 6 7, 6 4, 3 4, 3 7))", + 0.01, + "LINESTRING (9 9, 6 7)"); + } + + public void testLinesPolygon2() + { + String a = "MULTILINESTRING ((2 3, 2 7), (9 1, 9 8, 4 9))"; + String b = "POLYGON ((3 7, 6 8, 8 2, 3 4, 3 7))"; + checkDistance(a, b, 0.01, "LINESTRING (9 8, 6.3 7.1)"); + checkHD(a, b, 0.01, "LINESTRING (2 3, 5.5 3)"); + } + + public void testPolygonLineInteriorPoint() + { + checkDistanceStartPtLen("POLYGON ((2 8, 8 2, 2 1, 2 8))", + "LINESTRING (6 5, 4 7, 0 0, 9 1)", + 0.001, + "LINESTRING (4.557 2.9937, 2.4 4.2)", 0.01); + } + + public void testPolygonPolygon() + { + String a = "POLYGON ((2 18, 18 18, 17 3, 2 2, 2 18))"; + String b = "POLYGON ((1 19, 5 12, 5 3, 14 10, 11 19, 19 19, 20 0, 1 1, 1 19))"; + checkDistance(b, a, 0.01, "LINESTRING (20 0, 17 3)"); + checkDistance(a, b, 0.01, "LINESTRING (6.6796875 18, 11 19)"); + checkHD(a, b, 0.01, "LINESTRING (6.6796875 18, 11 19)"); + } + + public void testPolygonPolygonHolesNested() + { + // B is contained in A + String a = "POLYGON ((1 19, 19 19, 19 1, 1 1, 1 19), (6 8, 11 14, 15 7, 6 8))"; + String b = "POLYGON ((2 18, 18 18, 18 2, 2 2, 2 18), (10 17, 3 7, 17 5, 10 17))"; + checkDistance(a, b, 0.01, "LINESTRING (9.8134765625 12.576171875, 7.8603515625 13.943359375)"); + checkDistance(b, a, 0.01, 0.0); + } + + public void testMultiPolygons() + { + String a = "MULTIPOLYGON (((1 1, 1 10, 5 1, 1 1)), ((4 17, 9 15, 9 6, 4 17)))"; + String b = "MULTIPOLYGON (((1 12, 4 13, 8 10, 1 12)), ((3 8, 7 7, 6 2, 3 8)))"; + checkDistance(a, b, 0.01, "LINESTRING (1 1, 5.4 3.2)"); + checkDistance(b, a, 0.001, "LINESTRING (2.669921875 12.556640625, 5.446115154109589 13.818546660958905)"); + } + + // Tests that target area interiors have distance = 0 + public void testLinePolygonCrossing() throws Exception + { + String wkt1 = "LINESTRING (2 5, 5 10, 6 4)"; + String wkt2 = "POLYGON ((1 9, 9 9, 9 1, 1 1, 1 9))"; + checkDistance(wkt1, wkt2, 0.01, "LINESTRING (5 10, 5 9)"); + } + + public void testNonVertexResult() + { + String wkt1 = "LINESTRING (1 1, 5 10, 9 1)"; + String wkt2 = "LINESTRING (0 10, 0 0, 10 0)"; + + checkHD(wkt1, wkt2, 0.01, "LINESTRING (6.5390625 6.537109375, 6.5390625 0)"); + checkDistance(wkt1, wkt2, 0.01, "LINESTRING (6.5390625 6.537109375, 6.5390625 0)"); + } + + public void testDirectedLines() throws Exception + { + String wkt1 = "LINESTRING (1 6, 3 5, 1 4)"; + String wkt2 = "LINESTRING (1 10, 9 5, 1 2)"; + checkDistance(wkt1, wkt2, 0.01, "LINESTRING (1 6, 2.797752808988764 8.876404494382022)"); + checkDistance(wkt2, wkt1, 0.01, "LINESTRING (9 5, 3 5)"); + } + + public void testDirectedLines2() throws Exception + { + String wkt1 = "LINESTRING (1 6, 3 5, 1 4)"; + String wkt2 = "LINESTRING (1 3, 1 9, 9 5, 1 1)"; + checkDistance(wkt1, wkt2, 0.01, "LINESTRING (3 5, 1 5)"); + checkDistance(wkt2, wkt1, 0.01, "LINESTRING (9 5, 3 5)"); + } + + //====================================================================== + + private static final double TOLERANCE = 0.001; + + private void checkHD(String wkt1, String wkt2, double tolerance, String wktExpected) + { + Geometry g1 = read(wkt1); + Geometry g2 = read(wkt2); + + Geometry result = DirectedHausdorffDistance.hausdorffDistanceLine(g1, g2, tolerance); + Geometry expected = read(wktExpected); + checkEqualExact(expected, result, TOLERANCE); + + /* + double resultDistance = DiscreteHausdorffDistance.distance(g1, g2); + double expectedDistance = expected.getLength(); + assertEquals(expectedDistance, resultDistance, TOLERANCE); + */ + } + + private void checkDistance(String wkt1, String wkt2, double tolerance, String wktExpected) { + Geometry g1 = read(wkt1); + Geometry g2 = read(wkt2); + + Geometry result = DirectedHausdorffDistance.distanceLine(g1, g2, tolerance); + Geometry expected = read(wktExpected); + checkEqualExact(expected, result, 100 * tolerance); + } + + private void checkDistanceStartPtLen(String wkt1, String wkt2, double tolerance, + String wktExpected, double resultTolerance) { + Geometry g1 = read(wkt1); + Geometry g2 = read(wkt2); + + Geometry result = DirectedHausdorffDistance.distanceLine(g1, g2, tolerance); + Geometry expected = read(wktExpected); + + Coordinate resultPt = result.getCoordinates()[0]; + Coordinate expectedPt = expected.getCoordinates()[0]; + checkEqualXY(expectedPt, resultPt, resultTolerance); + + double distResult = result.getLength(); + double distExpected = expected.getLength(); + assertEquals(distExpected, distResult, resultTolerance); + } + + private void checkDistance(String wkt1, String wkt2, double tolerance, + double expectedDistance) { + Geometry g1 = read(wkt1); + Geometry g2 = read(wkt2); + + Geometry result = DirectedHausdorffDistance.distanceLine(g1, g2, tolerance); + + double distResult = result.getLength(); + assertEquals(expectedDistance, distResult, TOLERANCE); + } + +} diff --git a/modules/core/src/test/java/test/jts/perf/operation/distance/HausdorffDistanceStressTest.java b/modules/core/src/test/java/test/jts/perf/operation/distance/HausdorffDistanceStressTest.java new file mode 100644 index 0000000000..028e341aa6 --- /dev/null +++ b/modules/core/src/test/java/test/jts/perf/operation/distance/HausdorffDistanceStressTest.java @@ -0,0 +1,98 @@ +/* + * Copyright (c) 2026 Martin Davis. + * + * All rights reserved. This program and the accompanying materials + * are made available under the terms of the Eclipse Public License 2.0 + * and Eclipse Distribution License v. 1.0 which accompanies this distribution. + * The Eclipse Public License is available at http://www.eclipse.org/legal/epl-v20.html + * and the Eclipse Distribution License is available at + * + * http://www.eclipse.org/org/documents/edl-v10.php. + */ +package test.jts.perf.operation.distance; + +import java.util.ArrayList; +import java.util.List; + +import org.locationtech.jts.algorithm.distance.DiscreteHausdorffDistance; +import org.locationtech.jts.algorithm.distance.DirectedHausdorffDistance; +import org.locationtech.jts.geom.Coordinate; +import org.locationtech.jts.geom.Geometry; +import org.locationtech.jts.geom.GeometryFactory; +import org.locationtech.jts.geom.Point; +import org.locationtech.jts.geom.util.SineStarFactory; +import org.locationtech.jts.util.Debug; + +public class HausdorffDistanceStressTest { + static final int MAX_ITER = 10000; + + static final boolean VERBOSE = false; + + static GeometryFactory geomFact = new GeometryFactory(); + + public static void main(String[] args) { + + HausdorffDistanceStressTest test = new HausdorffDistanceStressTest(); + test.run(); + } + + private void run() { + Geometry b = createSineStar(1000, 300, 50); + Debug.println(b); + + List grid = createCircleGrid(2000, 100); + + int iter = 0; + for (Geometry g : grid) { + iter++; + Debug.println(iter + " ----------------"); + Debug.println(g); + checkHausdorff(g, b); + + //if (iter > 10) break; + } + } + + private void checkHausdorff(Geometry g, Geometry b) { + double distDHD = DiscreteHausdorffDistance.distance(g, b, 0.01); + + double distHD = DirectedHausdorffDistance.hausdorffDistanceLine(g, b, 0.01).getLength(); + //-- performance testing only + //double distDHD = distHD; + + double err = Math.abs(distDHD - distHD) / (distDHD + distHD); + + Debug.println(distDHD + " " + distHD + " err = " + err); + if (err > .01) { + System.out.println("<<<<<<<<<<< ERROR!"); + } + } + + private List createCircleGrid(double size, int nSide) { + List geoms = new ArrayList(); + + double inc = size / nSide; + + for (int i = 0; i < nSide; i++) { + for (int j = 0; j < nSide; j++) { + Coordinate p = new Coordinate(i * inc, j * inc); + Point pt = geomFact.createPoint(p); + Geometry buf = pt.buffer(inc); + geoms.add(buf); + } + } + return geoms; + } + + Geometry createSineStar(double loc, double size, int nPts) + { + SineStarFactory gsf = new SineStarFactory(geomFact); + gsf.setCentre(new Coordinate(loc, loc)); + gsf.setSize(size); + gsf.setNumPoints(nPts); + + Geometry g = gsf.createSineStar().getBoundary(); + + return g; + } +} From 5fc963adaaa1b47c1213712bdec852c4aa51ed4b Mon Sep 17 00:00:00 2001 From: Martin Davis Date: Fri, 30 Jan 2026 14:17:55 -0800 Subject: [PATCH 02/43] Expose IndexedPointInPolygonsLocater --- .../jts/algorithm/construct/IndexedDistanceToPoint.java | 1 + .../jts/algorithm/distance/DirectedHausdorffDistance.java | 2 +- .../IndexedPointInPolygonsLocator.java | 7 ++++--- 3 files changed, 6 insertions(+), 4 deletions(-) rename modules/core/src/main/java/org/locationtech/jts/algorithm/{construct => locate}/IndexedPointInPolygonsLocator.java (90%) diff --git a/modules/core/src/main/java/org/locationtech/jts/algorithm/construct/IndexedDistanceToPoint.java b/modules/core/src/main/java/org/locationtech/jts/algorithm/construct/IndexedDistanceToPoint.java index cac5810db9..9ba6f208dc 100644 --- a/modules/core/src/main/java/org/locationtech/jts/algorithm/construct/IndexedDistanceToPoint.java +++ b/modules/core/src/main/java/org/locationtech/jts/algorithm/construct/IndexedDistanceToPoint.java @@ -11,6 +11,7 @@ */ package org.locationtech.jts.algorithm.construct; +import org.locationtech.jts.algorithm.locate.IndexedPointInPolygonsLocator; import org.locationtech.jts.geom.Coordinate; import org.locationtech.jts.geom.Geometry; import org.locationtech.jts.geom.Location; diff --git a/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java b/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java index a6d25a079b..fc628e9671 100644 --- a/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java +++ b/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java @@ -14,8 +14,8 @@ import java.util.Iterator; import java.util.PriorityQueue; -import org.locationtech.jts.algorithm.construct.IndexedPointInPolygonsLocator; import org.locationtech.jts.algorithm.construct.LargestEmptyCircle; +import org.locationtech.jts.algorithm.locate.IndexedPointInPolygonsLocator; import org.locationtech.jts.geom.Coordinate; import org.locationtech.jts.geom.Dimension; import org.locationtech.jts.geom.Geometry; diff --git a/modules/core/src/main/java/org/locationtech/jts/algorithm/construct/IndexedPointInPolygonsLocator.java b/modules/core/src/main/java/org/locationtech/jts/algorithm/locate/IndexedPointInPolygonsLocator.java similarity index 90% rename from modules/core/src/main/java/org/locationtech/jts/algorithm/construct/IndexedPointInPolygonsLocator.java rename to modules/core/src/main/java/org/locationtech/jts/algorithm/locate/IndexedPointInPolygonsLocator.java index a5745c5102..68ec189846 100644 --- a/modules/core/src/main/java/org/locationtech/jts/algorithm/construct/IndexedPointInPolygonsLocator.java +++ b/modules/core/src/main/java/org/locationtech/jts/algorithm/locate/IndexedPointInPolygonsLocator.java @@ -9,12 +9,10 @@ * * http://www.eclipse.org/org/documents/edl-v10.php. */ -package org.locationtech.jts.algorithm.construct; +package org.locationtech.jts.algorithm.locate; import java.util.List; -import org.locationtech.jts.algorithm.locate.IndexedPointInAreaLocator; -import org.locationtech.jts.algorithm.locate.PointOnGeometryLocator; import org.locationtech.jts.geom.Coordinate; import org.locationtech.jts.geom.Envelope; import org.locationtech.jts.geom.Geometry; @@ -24,9 +22,12 @@ /** * Determines the location of a point in the polygonal elements of a geometry. + * The polygons may overlap. * Uses spatial indexing to provide efficient performance. * * @author mdavis + * + * @see IndexedPointInAreaLocator * */ public class IndexedPointInPolygonsLocator implements PointOnGeometryLocator { From ee3d24de1a2f3fb440c0173da7cad3fb5179f049 Mon Sep 17 00:00:00 2001 From: Martin Davis Date: Fri, 30 Jan 2026 14:28:27 -0800 Subject: [PATCH 03/43] Javadoc --- .../jts/algorithm/distance/DirectedHausdorffDistance.java | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java b/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java index fc628e9671..380eba8f8e 100644 --- a/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java +++ b/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java @@ -51,7 +51,8 @@ *

* Points, lines and polygons are supported as input. * If the query geometry is polygonal, - * the point of maximum distance may occur in the interior of a polygon. + * the point at maximum distance may occur in the interior of a query polygon. + * For a polygonal target geometry the point always lies on the boundary. *

* The class can be used in prepared mode. * Creating an instance on a target geometry caches indexes for that geometry. From 6f438f541255227d14dc84c943f6861a7a10d011 Mon Sep 17 00:00:00 2001 From: Martin Davis Date: Fri, 30 Jan 2026 16:50:45 -0800 Subject: [PATCH 04/43] Add short circuiting for isFullyWithinDistance --- .../jtstest/function/DistanceFunctions.java | 6 +- .../jtstest/function/SelectionFunctions.java | 15 +++- .../distance/DirectedHausdorffDistance.java | 89 +++++++++++++------ .../DirectedHausdorffDistanceTest.java | 23 +++++ 4 files changed, 101 insertions(+), 32 deletions(-) diff --git a/modules/app/src/main/java/org/locationtech/jtstest/function/DistanceFunctions.java b/modules/app/src/main/java/org/locationtech/jtstest/function/DistanceFunctions.java index af99d66d36..2926529470 100644 --- a/modules/app/src/main/java/org/locationtech/jtstest/function/DistanceFunctions.java +++ b/modules/app/src/main/java/org/locationtech/jtstest/function/DistanceFunctions.java @@ -120,7 +120,7 @@ public static Geometry nearestPointsIndexedEachB(Geometry a, Geometry b) { } @Metadata(description="Directed Hausdorff distance from A to B, up to tolerance") - public static double ohd(Geometry a, Geometry b, + public static double dhd(Geometry a, Geometry b, @Metadata(title="Distance tolerance") double distTol) { @@ -128,14 +128,14 @@ public static double ohd(Geometry a, Geometry b, } @Metadata(description="Directed Hausdorff distance from A to B, up to tolerance") - public static Geometry ohdLine(Geometry a, Geometry b, + public static Geometry dhdLine(Geometry a, Geometry b, @Metadata(title="Distance tolerance") double distTol) { return DirectedHausdorffDistance.distanceLine(a, b, distTol); } - @Metadata(description="Hausdorff distance from A to B, up to tolerance") + @Metadata(description="Hausdorff distance between A and B, up to tolerance") public static Geometry hdLine(Geometry a, Geometry b, @Metadata(title="Distance tolerance") double distTol) diff --git a/modules/app/src/main/java/org/locationtech/jtstest/function/SelectionFunctions.java b/modules/app/src/main/java/org/locationtech/jtstest/function/SelectionFunctions.java index c24f5d136a..cc212f9ca9 100644 --- a/modules/app/src/main/java/org/locationtech/jtstest/function/SelectionFunctions.java +++ b/modules/app/src/main/java/org/locationtech/jtstest/function/SelectionFunctions.java @@ -243,8 +243,19 @@ public static Geometry fullyWithinDistance(Geometry a, final Geometry mask, doub { return select(a, new GeometryPredicate() { public boolean isTrue(Geometry g) { - double tolerance = maximumDistance / 1e5; - return DirectedHausdorffDistance.isFullyWithinDistance(a, mask, maximumDistance, tolerance); + double tolerance = maximumDistance / 1e4; + return DirectedHausdorffDistance.isFullyWithinDistance(g, mask, maximumDistance, tolerance); + } + }); + } + + public static Geometry fullyWithinDistancePrep(Geometry a, final Geometry mask, double maximumDistance) + { + DirectedHausdorffDistance dhd = new DirectedHausdorffDistance(mask); + return select(a, new GeometryPredicate() { + public boolean isTrue(Geometry g) { + double tolerance = maximumDistance / 1e4; + return dhd.isFullyWithinDistance(g, maximumDistance, tolerance); } }); } diff --git a/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java b/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java index 380eba8f8e..46da092a8b 100644 --- a/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java +++ b/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java @@ -54,23 +54,25 @@ * the point at maximum distance may occur in the interior of a query polygon. * For a polygonal target geometry the point always lies on the boundary. *

+ * A common use case is to test whether a geometry A lies fully within a given + * distance of another one B. + * {@link #isFullyWithinDistance(Geometry, double, double)} + * can be used to test this efficiently. + * It implements heuristic checks and short-circuiting to improve performance. + * This can much more efficient than computing whether A is covered by a buffer of B. + *

* The class can be used in prepared mode. * Creating an instance on a target geometry caches indexes for that geometry. * Then {@link #computeDistancePoints(Geometry, double) * or {@link #isFullyWithinDistance(Geometry, double, double)} * can be called efficiently for multiple query geometries. *

- * A use case is to test whether a geometry A lies fully within a given - * distance of another one B. - * Using {@link #isFullyWithinDistance(Geometry, double, double)} - * is much more efficient than computing whether A is covered by a buffer of B. - *

* Due to the nature of the Hausdorff distance, * performance is not very sensitive to the distance tolerance, * so using a small tolerance is recommended. *

* This algorithm is easier to use, more accurate, - * and much faster than {@link DiscreteHausdorrffDistance}. + * and much faster than {@link DiscreteHausdorffDistance}. * * @author Martin Davis * @@ -113,6 +115,18 @@ public static LineString hausdorffDistanceLine(Geometry a, Geometry b, double to return a.getFactory().createLineString(pts); } + private static double distance(Coordinate[] pts) { + return pts[0].distance(pts[1]); + } + + private static Coordinate[] copyPair(Coordinate p0, Coordinate p1) { + return new Coordinate[] { p0.copy(), p1.copy() }; + } + + private static Coordinate[] copyReverse(Coordinate[] pts) { + return new Coordinate[] { pts[1].copy(), pts[0].copy() }; + } + private Geometry geomB; private IndexedFacetDistance distToB; private boolean isAreaB; @@ -128,21 +142,25 @@ public DirectedHausdorffDistance(Geometry b) { } public boolean isFullyWithinDistance(Geometry a, double maxDistance, double tolerance) { - //TODO: optimize with short-circuiting - Coordinate[] maxDistCoords = computeDistancePoints(a, tolerance); + //TODO: check envelopes + Coordinate[] maxDistCoords = computeDistancePoints(a, tolerance, maxDistance); return distance(maxDistCoords) <= maxDistance; } - - private static double distance(Coordinate[] pts) { - return pts[0].distance(pts[1]); - } private Coordinate[] computeDistancePoints(Geometry geomA, double tolerance) { + return computeDistancePoints(geomA, tolerance, -1.0); + } + + private Coordinate[] computeDistancePoints(Geometry geomA, double tolerance, double maxDistanceLimit) { if (geomA.getDimension() == Dimension.P) { - return computeAtPoints(geomA); + return computeAtPoints(geomA, maxDistanceLimit); } //TODO: handle mixed geoms with points - Coordinate[] maxDistPtsEdge = computeAtEdges(geomA, tolerance); + Coordinate[] maxDistPtsEdge = computeAtEdges(geomA, tolerance, maxDistanceLimit); + + if (isBeyondLimit(maxDistanceLimit, distance(maxDistPtsEdge))) { + return maxDistPtsEdge; + } /** * Polygonal query geometry may have an interior point as the farthest point. @@ -186,12 +204,14 @@ private Coordinate[] computeAtAreaInterior(Geometry geomA, double tolerance) { return null; } + //TODO: add short-circuiting based on maxDistanceLimit? + Point centerPt = LargestEmptyCircle.getCenter(geomB, polygonalA, tolerance); Coordinate ptA = centerPt.getCoordinate(); /** * If LEC centre is in B, the max distance is zero, so return null. * This will cause the computed segment distance to be returned, - * which is better since it occurs on a vertex (or edge?) + * which is preferred since it occurs on a vertex or edge. */ if (isInteriorB(ptA)) { return null; @@ -200,7 +220,7 @@ private Coordinate[] computeAtAreaInterior(Geometry geomA, double tolerance) { return copyReverse(nearPtsBA); } - private Coordinate[] computeAtPoints(Geometry geomA) { + private Coordinate[] computeAtPoints(Geometry geomA, double maxDistanceLimit) { double maxDist = -1.0;; Coordinate[] maxDistPtsAB = null; Iterator geomi = new GeometryCollectionIterator(geomA); @@ -223,19 +243,14 @@ private Coordinate[] computeAtPoints(Geometry geomA) { maxDist = dist; maxDistPtsAB = copyReverse(nearPtsBA); } + if (isBeyondLimit(maxDistanceLimit, maxDist)) { + break; + } } return maxDistPtsAB; } - - private static Coordinate[] copyPair(Coordinate p0, Coordinate p1) { - return new Coordinate[] { p0.copy(), p1.copy() }; - } - - private static Coordinate[] copyReverse(Coordinate[] pts) { - return new Coordinate[] { pts[1].copy(), pts[0].copy() }; - } - public Coordinate[] computeAtEdges(Geometry geomA, double tolerance) { + public Coordinate[] computeAtEdges(Geometry geomA, double tolerance, double maxDistanceLimit) { PriorityQueue segQueue = createSegQueue(geomA); DHDSegment ohdSeg = null; @@ -245,8 +260,10 @@ public Coordinate[] computeAtEdges(Geometry geomA, double tolerance) { iter++; // get the segment with greatest distance maxDistSeg = segQueue.remove(); - - if (maxDistSeg.length() <= tolerance) { + + double maxDist = maxDistSeg.maxDistance(); + if (maxDistSeg.length() <= tolerance + || isWithinLimit(maxDistanceLimit, maxDistSeg.maxDistanceBound())) { ohdSeg = maxDistSeg; break; } @@ -268,6 +285,14 @@ public Coordinate[] computeAtEdges(Geometry geomA, double tolerance) { return copyPair(maxPt, maxPt); } + private boolean isBeyondLimit(double maxDistanceLimit, double maxDist) { + return maxDistanceLimit >= 0 && maxDist > maxDistanceLimit; + } + + private boolean isWithinLimit(double maxDistanceLimit, double maxDist) { + return maxDistanceLimit >= 0 && maxDist <= maxDistanceLimit; + } + private void addNonInterior(DHDSegment ohdSegment, PriorityQueue segQueue) { //-- discard segment if it is interior to a polygon if (isInterior(ohdSegment)) { @@ -402,6 +427,16 @@ public double length() { return p0.distance(p1); } + public double maxDistance() { + double dist0 = p0.distance(nearPt0); + double dist1 = p1.distance(nearPt1); + return Math.max(dist0,dist1); + } + + public double maxDistanceBound() { + return maxDistanceBound; + } + public Coordinate[] getMaxDistPts() { double dist0 = p0.distance(nearPt0); double dist1 = p1.distance(nearPt1); diff --git a/modules/core/src/test/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistanceTest.java b/modules/core/src/test/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistanceTest.java index 5d17932b33..cbaf8a91f8 100644 --- a/modules/core/src/test/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistanceTest.java +++ b/modules/core/src/test/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistanceTest.java @@ -146,6 +146,18 @@ public void testDirectedLines2() throws Exception checkDistance(wkt2, wkt1, 0.01, "LINESTRING (9 5, 3 5)"); } + //----------------------------------------------------- + + public void testFullyWithinDistancePolygons() throws Exception + { + String a = "POLYGON ((1 4, 4 4, 4 1, 1 1, 1 4))"; + String b = "POLYGON ((10 10, 10 15, 15 15, 15 10, 10 10))"; + checkFullyWithinDistance(a, b, 5, 0.01, false); + checkFullyWithinDistance(a, b, 10, 0.01, false); + checkFullyWithinDistance(a, b, 20, 0.01, true); + } + + //====================================================================== private static final double TOLERANCE = 0.001; @@ -203,4 +215,15 @@ private void checkDistance(String wkt1, String wkt2, double tolerance, assertEquals(expectedDistance, distResult, TOLERANCE); } + + private void checkFullyWithinDistance(String a, String b, double distance, double tolerance, boolean expected) { + Geometry g1 = read(a); + Geometry g2 = read(b); + + boolean result = DirectedHausdorffDistance.isFullyWithinDistance(g1, g2, distance, tolerance); + + assertEquals(expected, result); + } + + } From b27bfddf529a65ae64cf48a8c51fc6c7a56912b4 Mon Sep 17 00:00:00 2001 From: Martin Davis Date: Sat, 31 Jan 2026 10:58:52 -0800 Subject: [PATCH 05/43] Add stress test for isFullyWithinDistance --- .../distance/HausdorffDistanceStressTest.java | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/modules/core/src/test/java/test/jts/perf/operation/distance/HausdorffDistanceStressTest.java b/modules/core/src/test/java/test/jts/perf/operation/distance/HausdorffDistanceStressTest.java index 028e341aa6..f7087c252d 100644 --- a/modules/core/src/test/java/test/jts/perf/operation/distance/HausdorffDistanceStressTest.java +++ b/modules/core/src/test/java/test/jts/perf/operation/distance/HausdorffDistanceStressTest.java @@ -66,6 +66,20 @@ private void checkHausdorff(Geometry g, Geometry b) { if (err > .01) { System.out.println("<<<<<<<<<<< ERROR!"); } + + checkFullyWithinDistance(g, b); + } + + private void checkFullyWithinDistance(Geometry g, Geometry b) { + double distDHD = DirectedHausdorffDistance.distanceLine(g, b, 0.01).getLength(); + double tol = distDHD / 1000; + boolean isWithin = DirectedHausdorffDistance.isFullyWithinDistance(g, b, 1.25 * distDHD, tol); + boolean isBeyond = ! DirectedHausdorffDistance.isFullyWithinDistance(g, b, 0.75 * distDHD, tol); + if (! (isWithin && isBeyond)) { + System.out.format("ioWithin = %b isBeyond = %b\n", isWithin, isBeyond); + System.out.println("<<<<<<<<<<< ERROR!"); + DirectedHausdorffDistance.isFullyWithinDistance(g, b, 0.75 * distDHD,tol); + } } private List createCircleGrid(double size, int nSide) { From 12bc572f54b136d868a237b62154183cb04f181d Mon Sep 17 00:00:00 2001 From: Martin Davis Date: Sat, 31 Jan 2026 11:22:11 -0800 Subject: [PATCH 06/43] Add early discard for interior segments --- .../distance/DirectedHausdorffDistance.java | 153 ++++++++++-------- .../DirectedHausdorffDistanceTest.java | 16 +- 2 files changed, 101 insertions(+), 68 deletions(-) diff --git a/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java b/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java index 46da092a8b..552d73e0a4 100644 --- a/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java +++ b/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java @@ -24,6 +24,7 @@ import org.locationtech.jts.geom.LineString; import org.locationtech.jts.geom.Location; import org.locationtech.jts.geom.Point; +import org.locationtech.jts.io.WKTWriter; import org.locationtech.jts.operation.distance.IndexedFacetDistance; /** @@ -158,7 +159,7 @@ private Coordinate[] computeDistancePoints(Geometry geomA, double tolerance, dou //TODO: handle mixed geoms with points Coordinate[] maxDistPtsEdge = computeAtEdges(geomA, tolerance, maxDistanceLimit); - if (isBeyondLimit(maxDistanceLimit, distance(maxDistPtsEdge))) { + if (isBeyondLimit(distance(maxDistPtsEdge), maxDistanceLimit)) { return maxDistPtsEdge; } @@ -175,51 +176,6 @@ && distance(maxDistPtsInterior) > distance(maxDistPtsEdge)) { return maxDistPtsEdge; } - /** - * If the query geometry A is polygonal, it is possible - * the farthest point lies in its interior. - * In this case it occurs at the centre of the Largest Empty Circle - * with B as obstacles and query geometry as constraint. - * - * This is a potentially expensive computation, - * especially for small tolerances. - * It is avoided where possible if heuristic checks - * can determine that the max distance is at - * the previously computed edge points. - * - * @param geomA - * @param tolerance - * @return the maximum distance point pair at an interior point of A, - * or null if it is known to not occur at an interior point - */ - private Coordinate[] computeAtAreaInterior(Geometry geomA, double tolerance) { - //TODO: extract polygonal geoms from A - Geometry polygonalA = geomA; - - /** - * Optimization - skip if A interior cannot intersect B, - * and thus farther point must lie on A segment - */ - if (polygonalA.getEnvelopeInternal().disjoint(geomB.getEnvelopeInternal())) { - return null; - } - - //TODO: add short-circuiting based on maxDistanceLimit? - - Point centerPt = LargestEmptyCircle.getCenter(geomB, polygonalA, tolerance); - Coordinate ptA = centerPt.getCoordinate(); - /** - * If LEC centre is in B, the max distance is zero, so return null. - * This will cause the computed segment distance to be returned, - * which is preferred since it occurs on a vertex or edge. - */ - if (isInteriorB(ptA)) { - return null; - } - Coordinate[] nearPtsBA = distToB.nearestPoints(ptA); - return copyReverse(nearPtsBA); - } - private Coordinate[] computeAtPoints(Geometry geomA, double maxDistanceLimit) { double maxDist = -1.0;; Coordinate[] maxDistPtsAB = null; @@ -243,17 +199,17 @@ private Coordinate[] computeAtPoints(Geometry geomA, double maxDistanceLimit) { maxDist = dist; maxDistPtsAB = copyReverse(nearPtsBA); } - if (isBeyondLimit(maxDistanceLimit, maxDist)) { + if (isBeyondLimit(maxDist, maxDistanceLimit)) { break; } } return maxDistPtsAB; } - public Coordinate[] computeAtEdges(Geometry geomA, double tolerance, double maxDistanceLimit) { + private Coordinate[] computeAtEdges(Geometry geomA, double tolerance, double maxDistanceLimit) { PriorityQueue segQueue = createSegQueue(geomA); - - DHDSegment ohdSeg = null; + + boolean maxDistSegFound = false; DHDSegment maxDistSeg = null; long iter = 0; while (! segQueue.isEmpty()) { @@ -261,10 +217,18 @@ public Coordinate[] computeAtEdges(Geometry geomA, double tolerance, double maxD // get the segment with greatest distance maxDistSeg = segQueue.remove(); - double maxDist = maxDistSeg.maxDistance(); + /** + * Exit if segment length is within tolerance, so distance is accurate enough. + * + * If maxDistanceLimit is specified, short-circuit if: + * - if segment distance bound is less than distance limit, no other segment can be farther + * - if a point of segment is farther than limit, isFulyWithin must be false + */ if (maxDistSeg.length() <= tolerance - || isWithinLimit(maxDistanceLimit, maxDistSeg.maxDistanceBound())) { - ohdSeg = maxDistSeg; + || isWithinLimit(maxDistSeg.maxDistanceBound(), maxDistanceLimit) + || isBeyondLimit(maxDistSeg.maxDistance(), maxDistanceLimit) + ) { + maxDistSegFound = true; break; } @@ -273,23 +237,24 @@ public Coordinate[] computeAtEdges(Geometry geomA, double tolerance, double maxD addNonInterior(bisects[0], segQueue); addNonInterior(bisects[1], segQueue); } - if (ohdSeg != null) - return ohdSeg.getMaxDistPts(); + if (maxDistSegFound) + return maxDistSeg.getMaxDistPts(); /** - * If no OHD segment was found, all were inside the target. + * No DHD segment was found. + * This must be because all were inside the target. * In this case distance is zero. * Return a single coordinate of the input as a representative point */ - Coordinate maxPt = maxDistSeg.getEndpoint(0); + Coordinate maxPt = geomA.getCoordinate(); return copyPair(maxPt, maxPt); } - - private boolean isBeyondLimit(double maxDistanceLimit, double maxDist) { + + private static boolean isBeyondLimit(double maxDist, double maxDistanceLimit) { return maxDistanceLimit >= 0 && maxDist > maxDistanceLimit; } - private boolean isWithinLimit(double maxDistanceLimit, double maxDist) { + private static boolean isWithinLimit(double maxDist, double maxDistanceLimit) { return maxDistanceLimit >= 0 && maxDist <= maxDistanceLimit; } @@ -325,6 +290,61 @@ private boolean isInteriorB(Coordinate p) { return ptInAreaB.locate(p) == Location.INTERIOR; } + /** + * If the query geometry A is polygonal, it is possible + * the farthest point lies in its interior. + * In this case it occurs at the centre of the Largest Empty Circle + * with B as obstacles and query geometry as constraint. + * + * This is a potentially expensive computation, + * especially for small tolerances. + * It is avoided where possible if heuristic checks + * can determine that the max distance is at + * the previously computed edge points. + * + * @param geomA + * @param tolerance + * @return the maximum distance point pair at an interior point of A, + * or null if it is known to not occur at an interior point + */ + private Coordinate[] computeAtAreaInterior(Geometry geomA, double tolerance) { + //TODO: extract polygonal geoms from A + Geometry polygonalA = geomA; + + /** + * Optimization - skip if A interior cannot intersect B, + * and thus farther point must lie on A segment + */ + if (polygonalA.getEnvelopeInternal().disjoint(geomB.getEnvelopeInternal())) { + return null; + } + + //TODO: add short-circuiting based on maxDistanceLimit? + + Point centerPt = LargestEmptyCircle.getCenter(geomB, polygonalA, tolerance); + Coordinate ptA = centerPt.getCoordinate(); + /** + * If LEC centre is in B, the max distance is zero, so return null. + * This will cause the computed segment distance to be returned, + * which is preferred since it occurs on a vertex or edge. + */ + if (isInteriorB(ptA)) { + return null; + } + Coordinate[] nearPtsBA = distToB.nearestPoints(ptA); + return copyReverse(nearPtsBA); + } + + /** + * Creates the priority queue for segments. + * Segments which are interior to a polygonal target geometry + * are not added to the queue. + * So if a query geometry is fully covered by the target + * the returned queue is empty. + * + * @param geomA + * @return the segment priority queue + */ private PriorityQueue createSegQueue(Geometry geomA) { PriorityQueue priq = new PriorityQueue(); geomA.apply(new GeometryComponentFilter() { @@ -332,7 +352,7 @@ private PriorityQueue createSegQueue(Geometry geomA) { @Override public void filter(Geometry geom) { if (geom instanceof LineString) { - initSegments(geom.getCoordinates(), priq); + addSegments(geom.getCoordinates(), priq); } } }); @@ -344,7 +364,7 @@ public void filter(Geometry geom) { * @param pts * @param priq */ - private void initSegments(Coordinate[] pts, PriorityQueue priq) { + private void addSegments(Coordinate[] pts, PriorityQueue priq) { DHDSegment prevSeg = null; for (int i = 0; i < pts.length - 1; i++) { DHDSegment seg; @@ -357,10 +377,9 @@ private void initSegments(Coordinate[] pts, PriorityQueue priq) { } prevSeg = seg; /** - * Delay interior check until splitting, - * to avoid computing more than needed. + * Don't add interior segments, since their distance must be zero. */ - priq.add(seg); + addNonInterior(seg, priq); //System.out.println(seg.distance()); } } @@ -481,5 +500,9 @@ public DHDSegment[] bisect(IndexedFacetDistance distToB) { public int compareTo(DHDSegment o) { return -Double.compare(maxDistanceBound, o.maxDistanceBound); } + + public String toString() { + return WKTWriter.toLineString(p0, p1); + } } } diff --git a/modules/core/src/test/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistanceTest.java b/modules/core/src/test/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistanceTest.java index cbaf8a91f8..5035e6becf 100644 --- a/modules/core/src/test/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistanceTest.java +++ b/modules/core/src/test/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistanceTest.java @@ -130,7 +130,7 @@ public void testNonVertexResult() checkDistance(wkt1, wkt2, 0.01, "LINESTRING (6.5390625 6.537109375, 6.5390625 0)"); } - public void testDirectedLines() throws Exception + public void testDirectedLines() { String wkt1 = "LINESTRING (1 6, 3 5, 1 4)"; String wkt2 = "LINESTRING (1 10, 9 5, 1 2)"; @@ -138,7 +138,7 @@ public void testDirectedLines() throws Exception checkDistance(wkt2, wkt1, 0.01, "LINESTRING (9 5, 3 5)"); } - public void testDirectedLines2() throws Exception + public void testDirectedLines2() { String wkt1 = "LINESTRING (1 6, 3 5, 1 4)"; String wkt2 = "LINESTRING (1 3, 1 9, 9 5, 1 1)"; @@ -146,9 +146,19 @@ public void testDirectedLines2() throws Exception checkDistance(wkt2, wkt1, 0.01, "LINESTRING (9 5, 3 5)"); } + /** + * Tests that segments are detected as interior even for a large tolerance. + */ + public void testInteriorSegmentsLargeTol() + { + String a = "POLYGON ((4 6, 5 6, 5 5, 4 5, 4 6))"; + String b = "POLYGON ((1 9, 9 9, 9 1, 1 1, 1 9))"; + checkDistance(a, b, 2.0, 0.0); + } + //----------------------------------------------------- - public void testFullyWithinDistancePolygons() throws Exception + public void testFullyWithinDistancePolygons() { String a = "POLYGON ((1 4, 4 4, 4 1, 1 1, 1 4))"; String b = "POLYGON ((10 10, 10 15, 15 15, 15 10, 10 10))"; From 1b14eeaa6586c75eca73ade237a23e9dd1bce101 Mon Sep 17 00:00:00 2001 From: Martin Davis Date: Sat, 31 Jan 2026 13:33:14 -0800 Subject: [PATCH 07/43] Rename methods --- .../distance/DirectedHausdorffDistance.java | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java b/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java index 552d73e0a4..d3c585cbac 100644 --- a/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java +++ b/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java @@ -111,7 +111,7 @@ public static LineString hausdorffDistanceLine(Geometry a, Geometry b, double to pts = ptsAB; } else { - pts = copyReverse(ptsBA); + pts = pairReverse(ptsBA); } return a.getFactory().createLineString(pts); } @@ -120,11 +120,11 @@ private static double distance(Coordinate[] pts) { return pts[0].distance(pts[1]); } - private static Coordinate[] copyPair(Coordinate p0, Coordinate p1) { + private static Coordinate[] pair(Coordinate p0, Coordinate p1) { return new Coordinate[] { p0.copy(), p1.copy() }; } - private static Coordinate[] copyReverse(Coordinate[] pts) { + private static Coordinate[] pairReverse(Coordinate[] pts) { return new Coordinate[] { pts[1].copy(), pts[0].copy() }; } @@ -193,11 +193,11 @@ private Coordinate[] computeAtPoints(Geometry geomA, double maxDistanceLimit) { //-- check for interior point if (isInterior) { dist = 0; - nearPtsBA = copyPair( pA, pA ); + nearPtsBA = pair( pA, pA ); } if (dist > maxDist) { maxDist = dist; - maxDistPtsAB = copyReverse(nearPtsBA); + maxDistPtsAB = pairReverse(nearPtsBA); } if (isBeyondLimit(maxDist, maxDistanceLimit)) { break; @@ -247,7 +247,7 @@ private Coordinate[] computeAtEdges(Geometry geomA, double tolerance, double max * Return a single coordinate of the input as a representative point */ Coordinate maxPt = geomA.getCoordinate(); - return copyPair(maxPt, maxPt); + return pair(maxPt, maxPt); } private static boolean isBeyondLimit(double maxDist, double maxDistanceLimit) { @@ -332,7 +332,7 @@ private Coordinate[] computeAtAreaInterior(Geometry geomA, double tolerance) { return null; } Coordinate[] nearPtsBA = distToB.nearestPoints(ptA); - return copyReverse(nearPtsBA); + return pairReverse(nearPtsBA); } /** @@ -460,9 +460,9 @@ public Coordinate[] getMaxDistPts() { double dist0 = p0.distance(nearPt0); double dist1 = p1.distance(nearPt1); if (dist0 > dist1) { - return copyPair(p0, nearPt0); + return pair(p0, nearPt0); } - return copyPair(p1, nearPt1); + return pair(p1, nearPt1); } /** From 1bf2555651bfb9b228268111dd1014778ac1cc14 Mon Sep 17 00:00:00 2001 From: Martin Davis Date: Sat, 31 Jan 2026 16:32:38 -0800 Subject: [PATCH 08/43] Add isFullyWithinDistance envelope checks --- .../distance/DirectedHausdorffDistance.java | 27 ++++++++++++++++++- .../DirectedHausdorffDistanceTest.java | 17 ++++++++++++ 2 files changed, 43 insertions(+), 1 deletion(-) diff --git a/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java b/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java index d3c585cbac..c7673b7a08 100644 --- a/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java +++ b/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java @@ -18,6 +18,7 @@ import org.locationtech.jts.algorithm.locate.IndexedPointInPolygonsLocator; import org.locationtech.jts.geom.Coordinate; import org.locationtech.jts.geom.Dimension; +import org.locationtech.jts.geom.Envelope; import org.locationtech.jts.geom.Geometry; import org.locationtech.jts.geom.GeometryCollectionIterator; import org.locationtech.jts.geom.GeometryComponentFilter; @@ -143,11 +144,35 @@ public DirectedHausdorffDistance(Geometry b) { } public boolean isFullyWithinDistance(Geometry a, double maxDistance, double tolerance) { - //TODO: check envelopes + //-- envelope check + if (isBeyond(a.getEnvelopeInternal(), geomB.getEnvelopeInternal(), maxDistance)) + return false; + Coordinate[] maxDistCoords = computeDistancePoints(a, tolerance, maxDistance); return distance(maxDistCoords) <= maxDistance; } + /** + * Tests if a geometry must have a point farther than the maximum distance + * using the geometry envelopes. + * + * @param envA + * @param envB + * @param maxDistance + * @return true if geometry A must have a far point from B + */ + private static boolean isBeyond(Envelope envA, Envelope envB, double maxDistance) { + /** + * At least one point of the geometry lies on each edge of the envelope, + * so if any edge is far from the closest edge of the B envelope + * there must be a point further than the max distance. + */ + return envA.getMinX() < envB.getMinX() - maxDistance + || envA.getMinY() < envB.getMinY() - maxDistance + || envA.getMaxX() > envB.getMaxX() + maxDistance + || envA.getMaxY() > envB.getMaxY() + maxDistance; + } + private Coordinate[] computeDistancePoints(Geometry geomA, double tolerance) { return computeDistancePoints(geomA, tolerance, -1.0); } diff --git a/modules/core/src/test/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistanceTest.java b/modules/core/src/test/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistanceTest.java index 5035e6becf..2a962b927a 100644 --- a/modules/core/src/test/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistanceTest.java +++ b/modules/core/src/test/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistanceTest.java @@ -158,6 +158,15 @@ public void testInteriorSegmentsLargeTol() //----------------------------------------------------- + public void testFullyWithinDistanceLines() + { + String a = "MULTILINESTRING ((1 1, 3 3), (7 7, 9 9))"; + String b = "MULTILINESTRING ((1 9, 1 5), (6 4, 8 2))"; + checkFullyWithinDistance(a, b, 1, 0.01, false); + checkFullyWithinDistance(a, b, 4, 0.01, false); + checkFullyWithinDistance(a, b, 6, 0.01, true); + } + public void testFullyWithinDistancePolygons() { String a = "POLYGON ((1 4, 4 4, 4 1, 1 1, 1 4))"; @@ -167,6 +176,14 @@ public void testFullyWithinDistancePolygons() checkFullyWithinDistance(a, b, 20, 0.01, true); } + public void testFullyWithinDistancePolygonsNestedWithHole() + { + String a = "POLYGON ((2 8, 8 8, 8 2, 2 2, 2 8))"; + String b = "POLYGON ((1 9, 9 9, 9 1, 1 1, 1 9), (3 7, 7 7, 7 3, 3 3, 3 7))"; + checkFullyWithinDistance(a, b, 1, 0.01, false); + checkFullyWithinDistance(a, b, 2, 0.01, true); + checkFullyWithinDistance(a, b, 3, 0.01, true); + } //====================================================================== From 3e28c1d9f653717cd86966d559391c1831e58339 Mon Sep 17 00:00:00 2001 From: Martin Davis Date: Sat, 31 Jan 2026 16:32:52 -0800 Subject: [PATCH 09/43] Refine stress test --- .../perf/operation/distance/HausdorffDistanceStressTest.java | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/modules/core/src/test/java/test/jts/perf/operation/distance/HausdorffDistanceStressTest.java b/modules/core/src/test/java/test/jts/perf/operation/distance/HausdorffDistanceStressTest.java index f7087c252d..a1513ffe9a 100644 --- a/modules/core/src/test/java/test/jts/perf/operation/distance/HausdorffDistanceStressTest.java +++ b/modules/core/src/test/java/test/jts/perf/operation/distance/HausdorffDistanceStressTest.java @@ -73,8 +73,8 @@ private void checkHausdorff(Geometry g, Geometry b) { private void checkFullyWithinDistance(Geometry g, Geometry b) { double distDHD = DirectedHausdorffDistance.distanceLine(g, b, 0.01).getLength(); double tol = distDHD / 1000; - boolean isWithin = DirectedHausdorffDistance.isFullyWithinDistance(g, b, 1.25 * distDHD, tol); - boolean isBeyond = ! DirectedHausdorffDistance.isFullyWithinDistance(g, b, 0.75 * distDHD, tol); + boolean isWithin = DirectedHausdorffDistance.isFullyWithinDistance(g, b, 1.05 * distDHD, tol); + boolean isBeyond = ! DirectedHausdorffDistance.isFullyWithinDistance(g, b, 0.95 * distDHD, tol); if (! (isWithin && isBeyond)) { System.out.format("ioWithin = %b isBeyond = %b\n", isWithin, isBeyond); System.out.println("<<<<<<<<<<< ERROR!"); From 57f6119f42fc53d6a173d150283c1f6fa5839a95 Mon Sep 17 00:00:00 2001 From: Martin Davis Date: Mon, 2 Feb 2026 15:55:49 -0800 Subject: [PATCH 10/43] Javadoc --- .../distance/DirectedHausdorffDistance.java | 94 ++++++++++++++++--- 1 file changed, 80 insertions(+), 14 deletions(-) diff --git a/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java b/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java index c7673b7a08..4c22ba2731 100644 --- a/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java +++ b/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java @@ -31,20 +31,21 @@ /** * Computes the directed Hausdorff distance from one geometry to another, * up to an approximation distance tolerance. - * The directed Hausdorff distance (DHD) is defined as: - * 

- * DHD(A,B) = maxa ∈ A (maxb ∈ B (distance(a, b) )
- *
- * The DHD is the maximum distance any point + * The directed Hausdorff distance is the maximum distance any point * on a query geometry A can be from a target geometry B. - * Equivalently, every point in the query geometry is within the DHD distance + * Equivalently, every point in the query geometry is within that distance * of the target geometry. * The class can compute a pair of points at which the DHD is obtained: * [ farthest A point, nearest B point ]. *

- * The DHD is asymmetric: DHD(A,B) may not be equal to DHD(B,A). + *The directed Hausdorff distance (DHD) is defined as: + * 

+ * DHD(A,B) = maxa ∈ A (maxb ∈ B (distance(a, b) )
+ *
+ *

+ * DHD is asymmetric: DHD(A,B) may not be equal to DHD(B,A). * Hence it is not a distance metric. - * The Hausdorff distance is is a symmetric distance metric: + * The Hausdorff distance is is a symmetric distance metric defined as: * 

  * HD(A,B) = max(DHD(A,B), DHD(B,A))
  *
@@ -61,7 +62,9 @@ * {@link #isFullyWithinDistance(Geometry, double, double)} * can be used to test this efficiently. * It implements heuristic checks and short-circuiting to improve performance. - * This can much more efficient than computing whether A is covered by a buffer of B. + * This can much more efficient than computing whether A is covered by B.buffer(distance). + * It is also more accurate, since constructed buffers + * are only linearized approximations to the true buffer. *

* The class can be used in prepared mode. * Creating an instance on a target geometry caches indexes for that geometry. @@ -76,29 +79,72 @@ * This algorithm is easier to use, more accurate, * and much faster than {@link DiscreteHausdorffDistance}. * + *

KNOWN ISSUES

+ * * @author Martin Davis * */ public class DirectedHausdorffDistance { + /** + * Computes the directed Hausdorff distance of a query geometry A from a target one B. + * + * @param a the query geometry + * @param b the target geometry + * @param tolerance the approximation distance tolerance + * @return the directed Hausdorff distance + */ public static double distance(Geometry a, Geometry b, double tolerance) { DirectedHausdorffDistance hd = new DirectedHausdorffDistance(b); return distance(hd.computeDistancePoints(a, tolerance)); } + /** + * Computes a pair of points which attain the directed Hausdorff distance + * of a query geometry A from a target one B. + * + * @param a the query geometry + * @param b the target geometry + * @param tolerance the approximation distance tolerance + * @return a pair of points [ptA, ptB] demonstrating the distance + */ public static LineString distanceLine(Geometry a, Geometry b, double tolerance) { DirectedHausdorffDistance hd = new DirectedHausdorffDistance(b); return a.getFactory().createLineString(hd.computeDistancePoints(a, tolerance)); } + /** + * Computes whether a query geometry lies fully within a give distance of a target geometry. + * Equivalently, detects whether any point of the query geometry is farther + * from the target than the specified distance. + * This is the case if DHD(A, B) > maxDistance. + * + * @param a the query geometry + * @param b the target geometry + * @param maxDistance the distance limit + * @param tolerance the approximation distance tolerance + * @return true if the query geometry lies fully within the distance of the target + */ public static boolean isFullyWithinDistance(Geometry a, Geometry b, double maxDistance, double tolerance) { DirectedHausdorffDistance hd = new DirectedHausdorffDistance(b); return hd.isFullyWithinDistance(a, maxDistance, tolerance); } + /** + * Computes a pair of points which attain the Hausdorff distance + * between two geometries. + * + * @param a a geometry + * @param b a geometry + * @param tolerance the approximation distance tolerance + * @return a pair of points [ptA, ptB] demonstrating the Hausdorff distance + */ public static LineString hausdorffDistanceLine(Geometry a, Geometry b, double tolerance) { DirectedHausdorffDistance hdAB = new DirectedHausdorffDistance(b); @@ -134,6 +180,11 @@ private static Coordinate[] pairReverse(Coordinate[] pts) { private boolean isAreaB; private IndexedPointInPolygonsLocator ptInAreaB; + /** + * Create a new instance for a target geometry B + * + * @param b the geometry to compute the distance from + */ public DirectedHausdorffDistance(Geometry b) { geomB = b; distToB = new IndexedFacetDistance(b); @@ -143,6 +194,17 @@ public DirectedHausdorffDistance(Geometry b) { } } + /** + * Computes whether a query geometry lies fully within a give distance of the target geometry. + * Equivalently, detects whether any point of the query geometry is farther + * from the target than the specified distance. + * This is the case if DHD(A, B) > maxDistance. + * + * @param a the query geometry + * @param maxDistance the distance limit + * @param tolerance the approximation distance tolerance + * @return true if the query geometry lies fully within the distance of the target + */ public boolean isFullyWithinDistance(Geometry a, double maxDistance, double tolerance) { //-- envelope check if (isBeyond(a.getEnvelopeInternal(), geomB.getEnvelopeInternal(), maxDistance)) @@ -173,7 +235,15 @@ private static boolean isBeyond(Envelope envA, Envelope envB, double maxDistance || envA.getMaxY() > envB.getMaxY() + maxDistance; } - private Coordinate[] computeDistancePoints(Geometry geomA, double tolerance) { + /** + * Computes a pair of points which attain the directed Hausdorff distance + * of a query geometry A from the target. + * + * @param geomA the query geometry + * @param tolerance the approximation distance tolerance + * @return a pair of points [ptA, ptB] demonstrating the distance + */ + public Coordinate[] computeDistancePoints(Geometry geomA, double tolerance) { return computeDistancePoints(geomA, tolerance, -1.0); } @@ -462,10 +532,6 @@ private static Coordinate nearestPt(Coordinate p, IndexedFacetDistance indexDist public Coordinate getEndpoint(int index) { return index == 0 ? p0 : p1; } - - public Coordinate getNearestPt(int index) { - return index == 0 ? nearPt0 : nearPt1; - } public double length() { return p0.distance(p1); From 9d62ddf2d363996245306d70553ae812beff414a Mon Sep 17 00:00:00 2001 From: Martin Davis Date: Tue, 3 Feb 2026 08:59:09 -0800 Subject: [PATCH 11/43] Check interior for segment endpoint distance --- .../distance/DirectedHausdorffDistance.java | 59 +++++++++++-------- .../DirectedHausdorffDistanceTest.java | 28 ++++++++- 2 files changed, 61 insertions(+), 26 deletions(-) diff --git a/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java b/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java index 4c22ba2731..d466119c70 100644 --- a/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java +++ b/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java @@ -226,8 +226,8 @@ public boolean isFullyWithinDistance(Geometry a, double maxDistance, double tole private static boolean isBeyond(Envelope envA, Envelope envB, double maxDistance) { /** * At least one point of the geometry lies on each edge of the envelope, - * so if any edge is far from the closest edge of the B envelope - * there must be a point further than the max distance. + * so if any edge is farther than maxDistance from the closest edge of the B envelope + * there must be a point at that distance (or further). */ return envA.getMinX() < envB.getMinX() - maxDistance || envA.getMinY() < envB.getMinY() - maxDistance @@ -328,7 +328,7 @@ private Coordinate[] computeAtEdges(Geometry geomA, double tolerance, double max } //-- not within tolerance, so bisect segment and keep searching - DHDSegment[] bisects = maxDistSeg.bisect(distToB); + DHDSegment[] bisects = maxDistSeg.bisect(distToB, ptInAreaB); addNonInterior(bisects[0], segQueue); addNonInterior(bisects[1], segQueue); } @@ -353,29 +353,33 @@ private static boolean isWithinLimit(double maxDist, double maxDistanceLimit) { return maxDistanceLimit >= 0 && maxDist <= maxDistanceLimit; } - private void addNonInterior(DHDSegment ohdSegment, PriorityQueue segQueue) { + private void addNonInterior(DHDSegment segment, PriorityQueue segQueue) { //-- discard segment if it is interior to a polygon - if (isInterior(ohdSegment)) { + if (isInterior(segment)) { return; } - segQueue.add(ohdSegment); + segQueue.add(segment); } /** * Tests if segment is fully in the interior of the target geometry polygons (if any). * - * @param ohdSegment + * @param segment * @return */ - private boolean isInterior(DHDSegment ohdSegment) { + private boolean isInterior(DHDSegment segment) { if (! isAreaB) return false; - double segDist = distToB.distance(ohdSegment.getEndpoint(0), ohdSegment.getEndpoint(1)); - //-- if segment touches B it is not in interior + if (segment.maxDistance() > 0.0) { + return false; + } + //-- compute distance to B linework + double segDist = distToB.distance(segment.getEndpoint(0), segment.getEndpoint(1)); + //-- if segment touches B linework it is not in interior if (segDist == 0) return false; //-- only need to test one point to check interior - Coordinate pt = ohdSegment.getEndpoint(0); + Coordinate pt = segment.getEndpoint(0); boolean isInterior = isInteriorB(pt); return isInterior; } @@ -464,11 +468,11 @@ private void addSegments(Coordinate[] pts, PriorityQueue priq) { for (int i = 0; i < pts.length - 1; i++) { DHDSegment seg; if (i == 0) { - seg = DHDSegment.create(pts[i], pts[i + 1], distToB); + seg = DHDSegment.create(pts[i], pts[i + 1], distToB, ptInAreaB); } else { //-- optimize by avoiding recomputing pt distance - seg = DHDSegment.create(prevSeg, pts[i + 1], distToB); + seg = DHDSegment.create(prevSeg, pts[i + 1], distToB, ptInAreaB); } prevSeg = seg; /** @@ -481,15 +485,15 @@ private void addSegments(Coordinate[] pts, PriorityQueue priq) { private static class DHDSegment implements Comparable { - public static DHDSegment create(Coordinate p0, Coordinate p1, IndexedFacetDistance indexDist) { + public static DHDSegment create(Coordinate p0, Coordinate p1, IndexedFacetDistance indexDist, IndexedPointInPolygonsLocator ptInAreaB) { DHDSegment seg = new DHDSegment(p0, p1); - seg.init(indexDist); + seg.init(indexDist, ptInAreaB); return seg; } - public static DHDSegment create(DHDSegment prevSeg, Coordinate p1, IndexedFacetDistance indexDist) { + public static DHDSegment create(DHDSegment prevSeg, Coordinate p1, IndexedFacetDistance indexDist, IndexedPointInPolygonsLocator ptInAreaB) { DHDSegment seg = new DHDSegment(prevSeg.p1, p1); - seg.init(prevSeg.nearPt1, indexDist); + seg.init(prevSeg.nearPt1, indexDist, ptInAreaB); return seg; } @@ -512,19 +516,24 @@ private DHDSegment(Coordinate p0, Coordinate nearPt0, Coordinate p1, Coordinate computeMaxDistanceBound(); } - private void init(IndexedFacetDistance indexDist) { - nearPt0 = nearestPt(p0, indexDist); - nearPt1 = nearestPt(p1, indexDist); + private void init(IndexedFacetDistance indexDist, IndexedPointInPolygonsLocator ptInAreaB) { + nearPt0 = nearestPt(p0, indexDist, ptInAreaB); + nearPt1 = nearestPt(p1, indexDist, ptInAreaB); computeMaxDistanceBound(); } - private void init(Coordinate nearest0, IndexedFacetDistance indexDist) { + private void init(Coordinate nearest0, IndexedFacetDistance indexDist, IndexedPointInPolygonsLocator ptInAreaB) { nearPt0 = nearest0; - nearPt1 = nearestPt(p1, indexDist); + nearPt1 = nearestPt(p1, indexDist, ptInAreaB); computeMaxDistanceBound(); } - private static Coordinate nearestPt(Coordinate p, IndexedFacetDistance indexDist) { + private static Coordinate nearestPt(Coordinate p, IndexedFacetDistance indexDist, IndexedPointInPolygonsLocator ptInAreaB) { + if (ptInAreaB != null) { + if (ptInAreaB.locate(p) != Location.EXTERIOR) { + return p; + } + } Coordinate[] nearestPts = indexDist.nearestPoints(p); return nearestPts[0]; } @@ -572,12 +581,12 @@ private void computeMaxDistanceBound() { maxDistanceBound = maxDist + length() / 2; } - public DHDSegment[] bisect(IndexedFacetDistance distToB) { + public DHDSegment[] bisect(IndexedFacetDistance distToB, IndexedPointInPolygonsLocator ptInAreaB) { Coordinate mid = new Coordinate( (p0.x + p1.x) / 2, (p0.y + p1.y) / 2 ); - Coordinate nearPtMid = nearestPt(mid, distToB); + Coordinate nearPtMid = nearestPt(mid, distToB, ptInAreaB); return new DHDSegment[] { new DHDSegment(p0, nearPt0, mid, nearPtMid ), new DHDSegment(mid, nearPtMid, p1, nearPt1) diff --git a/modules/core/src/test/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistanceTest.java b/modules/core/src/test/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistanceTest.java index 2a962b927a..5082e2cb2c 100644 --- a/modules/core/src/test/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistanceTest.java +++ b/modules/core/src/test/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistanceTest.java @@ -42,6 +42,13 @@ public void testMultiPoints() checkHD(a, b, 0.01, "LINESTRING (6 6, 5 0)"); } + public void testPointPolygonInterior() + { + checkDistance("POINT (3 4)", "POLYGON ((1 9, 9 9, 9 1, 1 1, 1 9))", + 0.01, + 0); + } + public void testLineSegments() { checkHD("LINESTRING (0 0, 2 0)", "LINESTRING (0 0, 2 1)", @@ -79,7 +86,15 @@ public void testLinesPolygon2() checkHD(a, b, 0.01, "LINESTRING (2 3, 5.5 3)"); } - public void testPolygonLineInteriorPoint() + public void testPolygonLineCrossingBoundaryResult() + { + checkDistance("POLYGON ((2 8, 8 2, 2 1, 2 8))", + "LINESTRING (6 5, 4 7, 0 0, 8 4)", + 0.001, + "LINESTRING (2 8, 3.9384615384615387 6.892307692307693)"); + } + + public void testPolygonLineCrossingInteriorPoint() { checkDistanceStartPtLen("POLYGON ((2 8, 8 2, 2 1, 2 8))", "LINESTRING (6 5, 4 7, 0 0, 9 1)", @@ -156,6 +171,17 @@ public void testInteriorSegmentsLargeTol() checkDistance(a, b, 2.0, 0.0); } + /** + * Tests that segment endpoint nearest points + * which are interior to B have distance 0 + */ + public void testInteriorSegmentsSameExterior() + { + String a = "POLYGON ((1 9, 3 9, 4 5, 5.05 9, 9 9, 9 1, 1 1, 1 9))"; + String b = "POLYGON ((1 9, 9 9, 9 1, 1 1, 1 9))"; + checkDistance(a, b, 1.0, 0.0); + } + //----------------------------------------------------- public void testFullyWithinDistanceLines() From 0c1bf913ccac3411c7331d8f9c0ecb850b125241 Mon Sep 17 00:00:00 2001 From: Martin Davis Date: Wed, 4 Feb 2026 07:46:02 -0800 Subject: [PATCH 12/43] Javadoc, tests --- .../distance/DirectedHausdorffDistance.java | 6 +- .../DirectedHausdorffDistanceTest.java | 29 ++++ .../distance/HausdorffDistancePerfTest.java | 138 ++++++++++++++++++ 3 files changed, 170 insertions(+), 3 deletions(-) create mode 100644 modules/core/src/test/java/test/jts/perf/operation/distance/HausdorffDistancePerfTest.java diff --git a/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java b/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java index d466119c70..88c61a6402 100644 --- a/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java +++ b/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java @@ -82,7 +82,7 @@ *

KNOWN ISSUES

*
    *
  • if the two geometries are identical or nearly so, - * performance is slower than desirable. + * performance may be slow. *
* @author Martin Davis * @@ -206,10 +206,10 @@ public DirectedHausdorffDistance(Geometry b) { * @return true if the query geometry lies fully within the distance of the target */ public boolean isFullyWithinDistance(Geometry a, double maxDistance, double tolerance) { - //-- envelope check + //-- envelope checks if (isBeyond(a.getEnvelopeInternal(), geomB.getEnvelopeInternal(), maxDistance)) return false; - + Coordinate[] maxDistCoords = computeDistancePoints(a, tolerance, maxDistance); return distance(maxDistCoords) <= maxDistance; } diff --git a/modules/core/src/test/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistanceTest.java b/modules/core/src/test/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistanceTest.java index 5082e2cb2c..228217b1ef 100644 --- a/modules/core/src/test/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistanceTest.java +++ b/modules/core/src/test/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistanceTest.java @@ -184,6 +184,35 @@ public void testInteriorSegmentsSameExterior() //----------------------------------------------------- + //-- shows withinDistance envelope check not triggering for disconnected A + public void testFullyWithinDistancePoints() + { + String a = "MULTIPOINT ((1 9), (9 1))"; + String b = "MULTIPOINT ((1 1), (9 9))"; + checkFullyWithinDistance(a, b, 1, 0.01, false); + checkFullyWithinDistance(a, b, 8.1, 0.01, true); + } + + public void testFullyWithinDistanceDisconnectedLines() + { + String a = "MULTILINESTRING ((1 9, 2 9), (8 1, 9 1))"; + String b = "LINESTRING (9 9, 1 1)"; + checkFullyWithinDistance(a, b, 1, 0.01, false); + checkFullyWithinDistance(a, b, 6, 0.01, true); + checkFullyWithinDistance(b, a, 1, 0.01, false); + checkFullyWithinDistance(b, a, 7.1, 0.01, true); + } + + public void testFullyWithinDistanceDisconnectedPolygons() + { + String a = "MULTIPOLYGON (((1 9, 2 9, 2 8, 1 8, 1 9)), ((8 2, 9 2, 9 1, 8 1, 8 2)))"; + String b = "POLYGON ((1 2, 9 9, 2 1, 1 2))"; + checkFullyWithinDistance(a, b, 1, 0.01, false); + checkFullyWithinDistance(a, b, 5.3, 0.01, true); + checkFullyWithinDistance(b, a, 1, 0.01, false); + checkFullyWithinDistance(b, a, 7.1, 0.01, true); + } + public void testFullyWithinDistanceLines() { String a = "MULTILINESTRING ((1 1, 3 3), (7 7, 9 9))"; diff --git a/modules/core/src/test/java/test/jts/perf/operation/distance/HausdorffDistancePerfTest.java b/modules/core/src/test/java/test/jts/perf/operation/distance/HausdorffDistancePerfTest.java new file mode 100644 index 0000000000..c6582afb13 --- /dev/null +++ b/modules/core/src/test/java/test/jts/perf/operation/distance/HausdorffDistancePerfTest.java @@ -0,0 +1,138 @@ +/* + * Copyright (c) 2026 Martin Davis. + * + * All rights reserved. This program and the accompanying materials + * are made available under the terms of the Eclipse Public License 2.0 + * and Eclipse Distribution License v. 1.0 which accompanies this distribution. + * The Eclipse Public License is available at http://www.eclipse.org/legal/epl-v20.html + * and the Eclipse Distribution License is available at + * + * http://www.eclipse.org/org/documents/edl-v10.php. + */ +package test.jts.perf.operation.distance; + +import java.util.ArrayList; +import java.util.List; + +import org.locationtech.jts.algorithm.distance.DiscreteHausdorffDistance; +import org.locationtech.jts.algorithm.distance.DirectedHausdorffDistance; +import org.locationtech.jts.geom.Coordinate; +import org.locationtech.jts.geom.Envelope; +import org.locationtech.jts.geom.Geometry; +import org.locationtech.jts.geom.GeometryFactory; +import org.locationtech.jts.geom.Point; +import org.locationtech.jts.geom.util.SineStarFactory; +import org.locationtech.jts.util.Debug; + +import test.jts.perf.PerformanceTestCase; +import test.jts.perf.PerformanceTestRunner; + +public class HausdorffDistancePerfTest extends PerformanceTestCase { + + static final int MAX_ITER = 10000; + + static final boolean VERBOSE = false; + + static GeometryFactory geomFact = new GeometryFactory(); + + public static void main(String[] args) { + PerformanceTestRunner.run(HausdorffDistancePerfTest.class); + } + + private Geometry b; + + private List grid; + + private int distance; + + private DirectedHausdorffDistance dhd; + + public HausdorffDistancePerfTest(String name) { + super(name); + setRunSize(new int[] { 100, 10_000, 40_000, 160_000, 1_000_000 }); + setRunIterations(1); + } + + public void startRun(int npts) + { + int nCircles = (int) Math.sqrt(npts); + System.out.println("\n------- Running with # circles = " + nCircles * nCircles); + b = createSineStar(1000, 300, 50); + dhd = new DirectedHausdorffDistance(b); + grid = createCircleGrid(2000, nCircles); + distance = 2000 / nCircles; + } + + public void runFullyWithin() { + Debug.println(b); + + int iter = 0; + for (Geometry g : grid) { + iter++; + Debug.println(iter + " ----------------"); + Debug.println(g); + //checkHausdorff(g, b); + //boolean isWithin = DirectedHausdorffDistance.isFullyWithinDistance(g, b, 4 * distance, 0.1); + boolean isWithin = dhd.isFullyWithinDistance(g, 4 * distance, 0.1); + + //if (iter > 10) break; + } + } + + private void checkHausdorff(Geometry g, Geometry b) { + double distDHD = DiscreteHausdorffDistance.distance(g, b, 0.01); + + double distHD = DirectedHausdorffDistance.hausdorffDistanceLine(g, b, 0.01).getLength(); + //-- performance testing only + //double distDHD = distHD; + + double err = Math.abs(distDHD - distHD) / (distDHD + distHD); + + Debug.println(distDHD + " " + distHD + " err = " + err); + if (err > .01) { + System.out.println("<<<<<<<<<<< ERROR!"); + } + + checkFullyWithinDistance(g, b); + } + + private void checkFullyWithinDistance(Geometry g, Geometry b) { + double distDHD = DirectedHausdorffDistance.distanceLine(g, b, 0.01).getLength(); + double tol = distDHD / 1000; + boolean isWithin = DirectedHausdorffDistance.isFullyWithinDistance(g, b, 1.05 * distDHD, tol); + boolean isBeyond = ! DirectedHausdorffDistance.isFullyWithinDistance(g, b, 0.95 * distDHD, tol); + if (! (isWithin && isBeyond)) { + System.out.format("ioWithin = %b isBeyond = %b\n", isWithin, isBeyond); + System.out.println("<<<<<<<<<<< ERROR!"); + DirectedHausdorffDistance.isFullyWithinDistance(g, b, 0.75 * distDHD,tol); + } + } + + private List createCircleGrid(double size, int nSide) { + List geoms = new ArrayList(); + + double inc = size / nSide; + + for (int i = 0; i < nSide; i++) { + for (int j = 0; j < nSide; j++) { + Coordinate p = new Coordinate(i * inc, j * inc); + Point pt = geomFact.createPoint(p); + Geometry buf = pt.buffer(inc); + geoms.add(buf); + } + } + return geoms; + } + + private Geometry createSineStar(double loc, double size, int nPts) + { + SineStarFactory gsf = new SineStarFactory(geomFact); + gsf.setCentre(new Coordinate(loc, loc)); + gsf.setSize(size); + gsf.setNumPoints(nPts); + + Geometry g = gsf.createSineStar().getBoundary(); + + return g; + } +} From dfacfd690c731148c7b3ddeb48ee49182ac55b37 Mon Sep 17 00:00:00 2001 From: Martin Davis Date: Thu, 5 Feb 2026 16:41:50 -0800 Subject: [PATCH 13/43] Factor out TargetDistance --- .../distance/DirectedHausdorffDistance.java | 119 ++++++++++-------- 1 file changed, 67 insertions(+), 52 deletions(-) diff --git a/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java b/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java index 88c61a6402..a8db347205 100644 --- a/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java +++ b/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java @@ -176,9 +176,7 @@ private static Coordinate[] pairReverse(Coordinate[] pts) { } private Geometry geomB; - private IndexedFacetDistance distToB; - private boolean isAreaB; - private IndexedPointInPolygonsLocator ptInAreaB; + private TargetDistance distanceToB; /** * Create a new instance for a target geometry B @@ -187,11 +185,7 @@ private static Coordinate[] pairReverse(Coordinate[] pts) { */ public DirectedHausdorffDistance(Geometry b) { geomB = b; - distToB = new IndexedFacetDistance(b); - isAreaB = b.getDimension() >= Dimension.A; - if (isAreaB) { - ptInAreaB = new IndexedPointInPolygonsLocator(b); - } + distanceToB = new TargetDistance(geomB); } /** @@ -281,10 +275,10 @@ private Coordinate[] computeAtPoints(Geometry geomA, double maxDistanceLimit) { continue; Coordinate pA = geomElemA.getCoordinate(); - Coordinate[] nearPtsBA = distToB.nearestPoints(pA); + Coordinate[] nearPtsBA = distanceToB.nearestPoints(pA); double dist = distance(nearPtsBA); - boolean isInterior = dist > 0 && isInteriorB(pA); + boolean isInterior = dist > 0 && distanceToB.isInterior(pA); //-- check for interior point if (isInterior) { dist = 0; @@ -328,7 +322,7 @@ private Coordinate[] computeAtEdges(Geometry geomA, double tolerance, double max } //-- not within tolerance, so bisect segment and keep searching - DHDSegment[] bisects = maxDistSeg.bisect(distToB, ptInAreaB); + DHDSegment[] bisects = maxDistSeg.bisect(distanceToB); addNonInterior(bisects[0], segQueue); addNonInterior(bisects[1], segQueue); } @@ -368,27 +362,12 @@ private void addNonInterior(DHDSegment segment, PriorityQueue segQue * @return */ private boolean isInterior(DHDSegment segment) { - if (! isAreaB) - return false; if (segment.maxDistance() > 0.0) { return false; } - //-- compute distance to B linework - double segDist = distToB.distance(segment.getEndpoint(0), segment.getEndpoint(1)); - //-- if segment touches B linework it is not in interior - if (segDist == 0) - return false; - //-- only need to test one point to check interior - Coordinate pt = segment.getEndpoint(0); - boolean isInterior = isInteriorB(pt); - return isInterior; + return distanceToB.isInterior(segment.getEndpoint(0), segment.getEndpoint(1)); } - private boolean isInteriorB(Coordinate p) { - if (! isAreaB) return false; - return ptInAreaB.locate(p) == Location.INTERIOR; - } - /** * If the query geometry A is polygonal, it is possible * the farthest point lies in its interior. @@ -427,10 +406,10 @@ private Coordinate[] computeAtAreaInterior(Geometry geomA, double tolerance) { * This will cause the computed segment distance to be returned, * which is preferred since it occurs on a vertex or edge. */ - if (isInteriorB(ptA)) { + if (distanceToB.isInterior(ptA)) { return null; } - Coordinate[] nearPtsBA = distToB.nearestPoints(ptA); + Coordinate[] nearPtsBA = distanceToB.nearestPoints(ptA); return pairReverse(nearPtsBA); } @@ -468,11 +447,11 @@ private void addSegments(Coordinate[] pts, PriorityQueue priq) { for (int i = 0; i < pts.length - 1; i++) { DHDSegment seg; if (i == 0) { - seg = DHDSegment.create(pts[i], pts[i + 1], distToB, ptInAreaB); + seg = DHDSegment.create(pts[i], pts[i + 1], distanceToB); } else { //-- optimize by avoiding recomputing pt distance - seg = DHDSegment.create(prevSeg, pts[i + 1], distToB, ptInAreaB); + seg = DHDSegment.create(prevSeg, pts[i + 1], distanceToB); } prevSeg = seg; /** @@ -483,17 +462,63 @@ private void addSegments(Coordinate[] pts, PriorityQueue priq) { } } + private static class TargetDistance { + private IndexedFacetDistance distance; + private boolean isArea; + private IndexedPointInPolygonsLocator ptInArea; + + public TargetDistance(Geometry geom) { + distance = new IndexedFacetDistance(geom); + isArea = geom.getDimension() >= Dimension.A; + if (isArea) { + ptInArea = new IndexedPointInPolygonsLocator(geom); + } + } + + public Coordinate[] nearestPoints(Coordinate p) { + return distance.nearestPoints(p); + } + + public boolean isInterior(Coordinate p) { + if (! isArea) return false; + return ptInArea.locate(p) == Location.INTERIOR; + } + + public boolean isInterior(Coordinate p0, Coordinate p1) { + if (! isArea) + return false; + //-- compute distance to B linework + double segDist = distance.distance(p0, p1); + //-- if segment touches B linework it is not in interior + if (segDist == 0) + return false; + //-- only need to test one point to check interior + boolean isInterior = isInterior(p0); + return isInterior; + } + + public Coordinate nearestPt(Coordinate p) { + if (ptInArea != null) { + if (ptInArea.locate(p) != Location.EXTERIOR) { + return p; + } + } + Coordinate[] nearestPts = distance.nearestPoints(p); + return nearestPts[0]; + } + } + private static class DHDSegment implements Comparable { - public static DHDSegment create(Coordinate p0, Coordinate p1, IndexedFacetDistance indexDist, IndexedPointInPolygonsLocator ptInAreaB) { + public static DHDSegment create(Coordinate p0, Coordinate p1, TargetDistance dist) { DHDSegment seg = new DHDSegment(p0, p1); - seg.init(indexDist, ptInAreaB); + seg.init(dist); return seg; } - public static DHDSegment create(DHDSegment prevSeg, Coordinate p1, IndexedFacetDistance indexDist, IndexedPointInPolygonsLocator ptInAreaB) { + public static DHDSegment create(DHDSegment prevSeg, Coordinate p1, TargetDistance dist) { DHDSegment seg = new DHDSegment(prevSeg.p1, p1); - seg.init(prevSeg.nearPt1, indexDist, ptInAreaB); + seg.init(prevSeg.nearPt1, dist); return seg; } @@ -516,28 +541,18 @@ private DHDSegment(Coordinate p0, Coordinate nearPt0, Coordinate p1, Coordinate computeMaxDistanceBound(); } - private void init(IndexedFacetDistance indexDist, IndexedPointInPolygonsLocator ptInAreaB) { - nearPt0 = nearestPt(p0, indexDist, ptInAreaB); - nearPt1 = nearestPt(p1, indexDist, ptInAreaB); + private void init(TargetDistance dist) { + nearPt0 = dist.nearestPt(p0); + nearPt1 = dist.nearestPt(p1); computeMaxDistanceBound(); } - private void init(Coordinate nearest0, IndexedFacetDistance indexDist, IndexedPointInPolygonsLocator ptInAreaB) { + private void init(Coordinate nearest0, TargetDistance dist) { nearPt0 = nearest0; - nearPt1 = nearestPt(p1, indexDist, ptInAreaB); + nearPt1 = dist.nearestPt(p1); computeMaxDistanceBound(); } - private static Coordinate nearestPt(Coordinate p, IndexedFacetDistance indexDist, IndexedPointInPolygonsLocator ptInAreaB) { - if (ptInAreaB != null) { - if (ptInAreaB.locate(p) != Location.EXTERIOR) { - return p; - } - } - Coordinate[] nearestPts = indexDist.nearestPoints(p); - return nearestPts[0]; - } - public Coordinate getEndpoint(int index) { return index == 0 ? p0 : p1; } @@ -581,12 +596,12 @@ private void computeMaxDistanceBound() { maxDistanceBound = maxDist + length() / 2; } - public DHDSegment[] bisect(IndexedFacetDistance distToB, IndexedPointInPolygonsLocator ptInAreaB) { + public DHDSegment[] bisect(TargetDistance dist) { Coordinate mid = new Coordinate( (p0.x + p1.x) / 2, (p0.y + p1.y) / 2 ); - Coordinate nearPtMid = nearestPt(mid, distToB, ptInAreaB); + Coordinate nearPtMid = dist.nearestPt(mid); return new DHDSegment[] { new DHDSegment(p0, nearPt0, mid, nearPtMid ), new DHDSegment(mid, nearPtMid, p1, nearPt1) From 45ba4729d2bbb5c5f601dc8d52f3c60ec1fc694e Mon Sep 17 00:00:00 2001 From: Martin Davis Date: Fri, 6 Feb 2026 13:53:59 -0800 Subject: [PATCH 14/43] Add unit tests --- .../distance/DirectedHausdorffDistanceTest.java | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/modules/core/src/test/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistanceTest.java b/modules/core/src/test/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistanceTest.java index 228217b1ef..58e660111d 100644 --- a/modules/core/src/test/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistanceTest.java +++ b/modules/core/src/test/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistanceTest.java @@ -26,14 +26,14 @@ public static void main(String args[]) { public DirectedHausdorffDistanceTest(String name) { super(name); } - public void testPoints() + public void testPointPoint() { checkHD("POINT (0 0)", "POINT (1 1)", 0.01, "LINESTRING (0 0, 1 1)"); } - public void testMultiPoints() + public void testPointsPoints() { String a = "MULTIPOINT ((0 1), (2 3), (4 5), (6 6))"; String b = "MULTIPOINT ((0.1 0), (1 0), (2 0), (3 0), (4 0), (5 0))"; @@ -49,6 +49,13 @@ public void testPointPolygonInterior() 0); } + public void testPointsPolygon() + { + checkDistance("MULTIPOINT ((4 3), (2 8), (8 5))", "POLYGON ((6 9, 6 4, 9 1, 1 1, 6 9))", + 0.01, + "LINESTRING (2 8, 4.426966292134832 6.48314606741573)"); + } + public void testLineSegments() { checkHD("LINESTRING (0 0, 2 0)", "LINESTRING (0 0, 2 1)", From 5bd5ecfe94fd7da97421e5bed6dc3a525d9c1440 Mon Sep 17 00:00:00 2001 From: Martin Davis Date: Sat, 7 Feb 2026 17:05:04 -0800 Subject: [PATCH 15/43] Simplify IndexedFacetDistance API --- .../distance/DirectedHausdorffDistance.java | 7 ++++--- .../distance/IndexedFacetDistance.java | 20 ++++++------------- 2 files changed, 10 insertions(+), 17 deletions(-) diff --git a/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java b/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java index a8db347205..60c0c1228b 100644 --- a/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java +++ b/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java @@ -476,7 +476,8 @@ public TargetDistance(Geometry geom) { } public Coordinate[] nearestPoints(Coordinate p) { - return distance.nearestPoints(p); + Coordinate pd = distance.nearestPoint(p); + return new Coordinate[] {pd, p}; } public boolean isInterior(Coordinate p) { @@ -503,8 +504,8 @@ public Coordinate nearestPt(Coordinate p) { return p; } } - Coordinate[] nearestPts = distance.nearestPoints(p); - return nearestPts[0]; + Coordinate nearestPt = distance.nearestPoint(p); + return nearestPt; } } diff --git a/modules/core/src/main/java/org/locationtech/jts/operation/distance/IndexedFacetDistance.java b/modules/core/src/main/java/org/locationtech/jts/operation/distance/IndexedFacetDistance.java index 4ee62e5dc9..5439f8ef01 100644 --- a/modules/core/src/main/java/org/locationtech/jts/operation/distance/IndexedFacetDistance.java +++ b/modules/core/src/main/java/org/locationtech/jts/operation/distance/IndexedFacetDistance.java @@ -162,35 +162,27 @@ public Coordinate[] nearestPoints(Geometry g) * @param p the point coordinate * @return the nearest point on the base geometry and the point */ - public Coordinate[] nearestPoints(Coordinate p) + public Coordinate nearestPoint(Coordinate p) { CoordinateSequence seq = new CoordinateArraySequence(new Coordinate[] { p }); FacetSequence fs2 = new FacetSequence(seq, 0); Object nearest = cachedTree.nearestNeighbour(fs2.getEnvelope(), fs2, FACET_SEQ_DIST); FacetSequence fs1 = (FacetSequence) nearest; - return fs1.nearestLocations(fs2); + return fs1.nearestLocations(fs2)[0]; } - public double distance(Coordinate p) { - return distance(nearestPoints(p)); + return p.distance(nearestPoint(p)); } - private static double distance(Coordinate[] pts) { - return pts[0].distance(pts[1]); - } - - public Coordinate[] nearestPoints(Coordinate p0, Coordinate p1) + public double distance(Coordinate p0, Coordinate p1) { CoordinateSequence seq = new CoordinateArraySequence(new Coordinate[] { p0, p1 }); FacetSequence fs2 = new FacetSequence(seq, 0, 2); Object nearest = cachedTree.nearestNeighbour(fs2.getEnvelope(), fs2, FACET_SEQ_DIST); FacetSequence fs1 = (FacetSequence) nearest; - return fs1.nearestLocations(fs2); - } - - public double distance(Coordinate p0, Coordinate p1) { - return distance(nearestPoints(p0, p1)); + Coordinate[] loc = fs1.nearestLocations(fs2); + return loc[0].distance(loc[1]); } /** From c4962b3b4dbde924dcbd75839c938a5831a19725 Mon Sep 17 00:00:00 2001 From: Martin Davis Date: Sun, 8 Feb 2026 10:33:34 -0800 Subject: [PATCH 16/43] Formatting --- .../distance/IndexedFacetDistance.java | 32 +++++++++---------- 1 file changed, 16 insertions(+), 16 deletions(-) diff --git a/modules/core/src/main/java/org/locationtech/jts/operation/distance/IndexedFacetDistance.java b/modules/core/src/main/java/org/locationtech/jts/operation/distance/IndexedFacetDistance.java index 5439f8ef01..c187ac474d 100644 --- a/modules/core/src/main/java/org/locationtech/jts/operation/distance/IndexedFacetDistance.java +++ b/modules/core/src/main/java/org/locationtech/jts/operation/distance/IndexedFacetDistance.java @@ -117,22 +117,6 @@ public IndexedFacetDistance(Geometry geom) { cachedTree = FacetSequenceTreeBuilder.build(geom); } - /** - * Computes the distance from the base geometry to - * the given geometry. - * - * @param g the geometry to compute the distance to - * - * @return the computed distance - */ - public double distance(Geometry g) - { - Object[] obj = nearestFacets(g); - FacetSequence fs1 = (FacetSequence) obj[0]; - FacetSequence fs2 = (FacetSequence) obj[1]; - return fs1.distance(fs2); - } - private Object[] nearestFacets(Geometry g) { STRtree tree2 = FacetSequenceTreeBuilder.build(g); Object[] obj = cachedTree.nearestNeighbour(tree2, @@ -171,6 +155,22 @@ public Coordinate nearestPoint(Coordinate p) return fs1.nearestLocations(fs2)[0]; } + /** + * Computes the distance from the base geometry to + * the given geometry. + * + * @param g the geometry to compute the distance to + * + * @return the computed distance + */ + public double distance(Geometry g) + { + Object[] obj = nearestFacets(g); + FacetSequence fs1 = (FacetSequence) obj[0]; + FacetSequence fs2 = (FacetSequence) obj[1]; + return fs1.distance(fs2); + } + public double distance(Coordinate p) { return p.distance(nearestPoint(p)); } From 1090f6f7c73b5bb2af21246783630fa8428b0ac6 Mon Sep 17 00:00:00 2001 From: Martin Davis Date: Sun, 8 Feb 2026 10:33:53 -0800 Subject: [PATCH 17/43] Simplify coordinate handling --- .../distance/DirectedHausdorffDistance.java | 69 ++++++++++--------- 1 file changed, 36 insertions(+), 33 deletions(-) diff --git a/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java b/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java index 60c0c1228b..008f0d925d 100644 --- a/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java +++ b/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java @@ -275,18 +275,18 @@ private Coordinate[] computeAtPoints(Geometry geomA, double maxDistanceLimit) { continue; Coordinate pA = geomElemA.getCoordinate(); - Coordinate[] nearPtsBA = distanceToB.nearestPoints(pA); - double dist = distance(nearPtsBA); + Coordinate pB = distanceToB.nearestPoint(pA); + double dist = pA.distance(pB); boolean isInterior = dist > 0 && distanceToB.isInterior(pA); //-- check for interior point if (isInterior) { dist = 0; - nearPtsBA = pair( pA, pA ); + pB = pA; } if (dist > maxDist) { maxDist = dist; - maxDistPtsAB = pairReverse(nearPtsBA); + maxDistPtsAB = pair(pA, pB); } if (isBeyondLimit(maxDist, maxDistanceLimit)) { break; @@ -299,12 +299,12 @@ private Coordinate[] computeAtEdges(Geometry geomA, double tolerance, double max PriorityQueue segQueue = createSegQueue(geomA); boolean maxDistSegFound = false; - DHDSegment maxDistSeg = null; + DHDSegment segMaxDist = null; long iter = 0; while (! segQueue.isEmpty()) { iter++; // get the segment with greatest distance - maxDistSeg = segQueue.remove(); + segMaxDist = segQueue.remove(); /** * Exit if segment length is within tolerance, so distance is accurate enough. @@ -313,21 +313,24 @@ private Coordinate[] computeAtEdges(Geometry geomA, double tolerance, double max * - if segment distance bound is less than distance limit, no other segment can be farther * - if a point of segment is farther than limit, isFulyWithin must be false */ - if (maxDistSeg.length() <= tolerance - || isWithinLimit(maxDistSeg.maxDistanceBound(), maxDistanceLimit) - || isBeyondLimit(maxDistSeg.maxDistance(), maxDistanceLimit) + if (segMaxDist.length() <= tolerance + || isWithinLimit(segMaxDist.maxDistanceBound(), maxDistanceLimit) + || isBeyondLimit(segMaxDist.maxDistance(), maxDistanceLimit) ) { maxDistSegFound = true; break; } + //System.out.println(segMaxDist); + + //-- not within tolerance, so bisect segment and keep searching - DHDSegment[] bisects = maxDistSeg.bisect(distanceToB); + DHDSegment[] bisects = segMaxDist.bisect(distanceToB); addNonInterior(bisects[0], segQueue); addNonInterior(bisects[1], segQueue); } if (maxDistSegFound) - return maxDistSeg.getMaxDistPts(); + return segMaxDist.getMaxDistPts(); /** * No DHD segment was found. @@ -409,8 +412,8 @@ private Coordinate[] computeAtAreaInterior(Geometry geomA, double tolerance) { if (distanceToB.isInterior(ptA)) { return null; } - Coordinate[] nearPtsBA = distanceToB.nearestPoints(ptA); - return pairReverse(nearPtsBA); + Coordinate ptB = distanceToB.nearestFacetPoint(ptA); + return pair(ptA, ptB); } /** @@ -463,21 +466,30 @@ private void addSegments(Coordinate[] pts, PriorityQueue priq) { } private static class TargetDistance { - private IndexedFacetDistance distance; + private IndexedFacetDistance distanceToFacets; private boolean isArea; private IndexedPointInPolygonsLocator ptInArea; public TargetDistance(Geometry geom) { - distance = new IndexedFacetDistance(geom); + distanceToFacets = new IndexedFacetDistance(geom); isArea = geom.getDimension() >= Dimension.A; if (isArea) { ptInArea = new IndexedPointInPolygonsLocator(geom); } } - public Coordinate[] nearestPoints(Coordinate p) { - Coordinate pd = distance.nearestPoint(p); - return new Coordinate[] {pd, p}; + public Coordinate nearestFacetPoint(Coordinate p) { + return distanceToFacets.nearestPoint(p); + } + + public Coordinate nearestPoint(Coordinate p) { + if (ptInArea != null) { + if (ptInArea.locate(p) != Location.EXTERIOR) { + return p; + } + } + Coordinate nearestPt = distanceToFacets.nearestPoint(p); + return nearestPt; } public boolean isInterior(Coordinate p) { @@ -489,7 +501,7 @@ public boolean isInterior(Coordinate p0, Coordinate p1) { if (! isArea) return false; //-- compute distance to B linework - double segDist = distance.distance(p0, p1); + double segDist = distanceToFacets.distance(p0, p1); //-- if segment touches B linework it is not in interior if (segDist == 0) return false; @@ -497,16 +509,7 @@ public boolean isInterior(Coordinate p0, Coordinate p1) { boolean isInterior = isInterior(p0); return isInterior; } - - public Coordinate nearestPt(Coordinate p) { - if (ptInArea != null) { - if (ptInArea.locate(p) != Location.EXTERIOR) { - return p; - } - } - Coordinate nearestPt = distance.nearestPoint(p); - return nearestPt; - } + } private static class DHDSegment implements Comparable { @@ -543,14 +546,14 @@ private DHDSegment(Coordinate p0, Coordinate nearPt0, Coordinate p1, Coordinate } private void init(TargetDistance dist) { - nearPt0 = dist.nearestPt(p0); - nearPt1 = dist.nearestPt(p1); + nearPt0 = dist.nearestPoint(p0); + nearPt1 = dist.nearestPoint(p1); computeMaxDistanceBound(); } private void init(Coordinate nearest0, TargetDistance dist) { nearPt0 = nearest0; - nearPt1 = dist.nearestPt(p1); + nearPt1 = dist.nearestPoint(p1); computeMaxDistanceBound(); } @@ -602,7 +605,7 @@ public DHDSegment[] bisect(TargetDistance dist) { (p0.x + p1.x) / 2, (p0.y + p1.y) / 2 ); - Coordinate nearPtMid = dist.nearestPt(mid); + Coordinate nearPtMid = dist.nearestPoint(mid); return new DHDSegment[] { new DHDSegment(p0, nearPt0, mid, nearPtMid ), new DHDSegment(mid, nearPtMid, p1, nearPt1) From 8e878e7a3cad26717be2baadbc37bf0e4a2aa677 Mon Sep 17 00:00:00 2001 From: Martin Davis Date: Mon, 9 Feb 2026 14:32:47 -0800 Subject: [PATCH 18/43] Add FacetLocation --- .../jts/operation/distance/FacetLocation.java | 40 ++++++++++++++ .../jts/operation/distance/FacetSequence.java | 52 +++++++++++++++++++ .../distance/IndexedFacetDistance.java | 35 ++++++++----- 3 files changed, 114 insertions(+), 13 deletions(-) create mode 100644 modules/core/src/main/java/org/locationtech/jts/operation/distance/FacetLocation.java diff --git a/modules/core/src/main/java/org/locationtech/jts/operation/distance/FacetLocation.java b/modules/core/src/main/java/org/locationtech/jts/operation/distance/FacetLocation.java new file mode 100644 index 0000000000..3b0d543b87 --- /dev/null +++ b/modules/core/src/main/java/org/locationtech/jts/operation/distance/FacetLocation.java @@ -0,0 +1,40 @@ +package org.locationtech.jts.operation.distance; + +import org.locationtech.jts.geom.Coordinate; +import org.locationtech.jts.geom.CoordinateSequence; +import org.locationtech.jts.geom.CoordinateSequences; + +public class FacetLocation { + private CoordinateSequence seq; + private int index; + private Coordinate pt; + + public FacetLocation(CoordinateSequence seq, int index, Coordinate pt) { + this.seq = seq; + this.index = index; + this.pt = pt; + } + + public boolean isSameSegment(FacetLocation f) { + return seq == f.seq && index == f.index; + } + + public int getIndex() { + return index; + } + + public Coordinate getEndPoint(int i) { + if (i == 0) + return seq.getCoordinate(index); + else + return seq.getCoordinate(index + 1); + } + + public int normalize(int index) { + if (index >= seq.size() - 1 + && CoordinateSequences.isRing(seq)) { + return 0; + } + return index; + } +} diff --git a/modules/core/src/main/java/org/locationtech/jts/operation/distance/FacetSequence.java b/modules/core/src/main/java/org/locationtech/jts/operation/distance/FacetSequence.java index 44c587189a..ce07747f55 100644 --- a/modules/core/src/main/java/org/locationtech/jts/operation/distance/FacetSequence.java +++ b/modules/core/src/main/java/org/locationtech/jts/operation/distance/FacetSequence.java @@ -15,6 +15,7 @@ import org.locationtech.jts.algorithm.Distance; import org.locationtech.jts.geom.Coordinate; import org.locationtech.jts.geom.CoordinateSequence; +import org.locationtech.jts.geom.CoordinateSequences; import org.locationtech.jts.geom.Envelope; import org.locationtech.jts.geom.LineSegment; @@ -151,6 +152,57 @@ else if (isPointOther) { return locs; } + public FacetLocation nearestLocation(Coordinate p) + { + boolean isPoint = isPoint(); + + if (isPoint) { + return new FacetLocation(pts, 0, pts.getCoordinate(0)); + } + + return nearestLocationOnLine(p); + } + + private FacetLocation nearestLocationOnLine(Coordinate pt) + { + double minDistance = Double.MAX_VALUE; + int index = -1; + Coordinate nearestPt = null; + + for (int i = start; i < end - 1; i++) { + Coordinate q0 = pts.getCoordinate(i); + Coordinate q1 = pts.getCoordinate(i + 1); + double dist = Distance.pointToSegment(pt, q0, q1); + if (dist < minDistance) { + minDistance = dist; + nearestPt = nearestPoint(pt, q0, q1); + index = i; + //-- segments are half-open, so final endpoint belongs to next segment + if (dist == 0.0 && pt.equals2D(q1)) { + index++; + //-- normalize index for a ring + index = normalize(pts, index); + } + if (minDistance <= 0.0) + break; + } + } + return new FacetLocation(pts, index, nearestPt); + } + + private static int normalize(CoordinateSequence pts, int index) { + if (index >= pts.size() - 1 + && CoordinateSequences.isRing(pts)) { + index = 0; + } + return index; + } + + private Coordinate nearestPoint(Coordinate pt, Coordinate q0, Coordinate q1) { + LineSegment seg = new LineSegment(q0, q1); + return seg.closestPoint(pt); + } + private double computeDistanceLineLine(FacetSequence facetSeq, Coordinate[] locs) { // both linear - compute minimum segment-segment distance diff --git a/modules/core/src/main/java/org/locationtech/jts/operation/distance/IndexedFacetDistance.java b/modules/core/src/main/java/org/locationtech/jts/operation/distance/IndexedFacetDistance.java index c187ac474d..3615ca6ae3 100644 --- a/modules/core/src/main/java/org/locationtech/jts/operation/distance/IndexedFacetDistance.java +++ b/modules/core/src/main/java/org/locationtech/jts/operation/distance/IndexedFacetDistance.java @@ -24,8 +24,8 @@ import org.locationtech.jts.index.strtree.STRtree; /** - * Computes the distance between the facets (segments and vertices) - * of two {@link Geometry}s + * Computes the distance to the facets (segments and vertices) + * of a {@link Geometry} * using a Branch-and-Bound algorithm. * The Branch-and-Bound algorithm operates over a * traversal of R-trees built @@ -125,11 +125,11 @@ private Object[] nearestFacets(Geometry g) { } /** - * Computes the nearest points on the base geometry + * Computes the nearest points on the target geometry * and another geometry. * * @param g the geometry to compute the nearest location to - * @return the nearest points on the base and the argument geometry + * @return the nearest points on the target and the argument geometry */ public Coordinate[] nearestPoints(Geometry g) { @@ -139,24 +139,33 @@ public Coordinate[] nearestPoints(Geometry g) return fs1.nearestLocations(fs2); } + public FacetLocation nearestLocation(Coordinate p) { + CoordinateSequence seq = new CoordinateArraySequence(new Coordinate[] { p }); + FacetSequence fs = new FacetSequence(seq, 0); + Object nearest = cachedTree.nearestNeighbour(fs.getEnvelope(), fs, FACET_SEQ_DIST); + FacetSequence fsN = (FacetSequence) nearest; + return fsN.nearestLocation(p); + + } + /** - * Computes the nearest points on the base geometry - * and a point. + * Computes the nearest point on the target geometry + * to a point. * * @param p the point coordinate - * @return the nearest point on the base geometry and the point + * @return the nearest point on the target geometry */ public Coordinate nearestPoint(Coordinate p) { CoordinateSequence seq = new CoordinateArraySequence(new Coordinate[] { p }); - FacetSequence fs2 = new FacetSequence(seq, 0); - Object nearest = cachedTree.nearestNeighbour(fs2.getEnvelope(), fs2, FACET_SEQ_DIST); - FacetSequence fs1 = (FacetSequence) nearest; - return fs1.nearestLocations(fs2)[0]; + FacetSequence fs = new FacetSequence(seq, 0); + Object nearest = cachedTree.nearestNeighbour(fs.getEnvelope(), fs, FACET_SEQ_DIST); + FacetSequence fsN = (FacetSequence) nearest; + return fsN.nearestLocations(fs)[0]; } /** - * Computes the distance from the base geometry to + * Computes the distance from the target geometry to * the given geometry. * * @param g the geometry to compute the distance to @@ -186,7 +195,7 @@ public double distance(Coordinate p0, Coordinate p1) } /** - * Tests whether the base geometry lies within + * Tests whether the target geometry lies within * a specified distance of the given geometry. * * @param g the geometry to test From 9057ed133952173925186de0a2066c18e16dae3c Mon Sep 17 00:00:00 2001 From: Martin Davis Date: Mon, 9 Feb 2026 14:33:05 -0800 Subject: [PATCH 19/43] Check for same segments --- .../distance/DirectedHausdorffDistance.java | 36 +++++++++++++++++-- 1 file changed, 34 insertions(+), 2 deletions(-) diff --git a/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java b/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java index 008f0d925d..8758c826b8 100644 --- a/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java +++ b/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java @@ -26,6 +26,7 @@ import org.locationtech.jts.geom.Location; import org.locationtech.jts.geom.Point; import org.locationtech.jts.io.WKTWriter; +import org.locationtech.jts.operation.distance.FacetLocation; import org.locationtech.jts.operation.distance.IndexedFacetDistance; /** @@ -320,10 +321,18 @@ private Coordinate[] computeAtEdges(Geometry geomA, double tolerance, double max maxDistSegFound = true; break; } - //System.out.println(segMaxDist); - + //-- check for equal or coincident segments + if (segMaxDist.maxDistance() == 0.0) { + if (isSameSegment(segMaxDist)) + continue; + //System.out.println(segMaxDist); + isSameSegment(segMaxDist); + + } + //System.out.println(segMaxDist); + //-- not within tolerance, so bisect segment and keep searching DHDSegment[] bisects = segMaxDist.bisect(distanceToB); addNonInterior(bisects[0], segQueue); @@ -342,6 +351,25 @@ private Coordinate[] computeAtEdges(Geometry geomA, double tolerance, double max return pair(maxPt, maxPt); } + private boolean isSameSegment(DHDSegment seg) { + FacetLocation f0 = distanceToB.nearestLocation(seg.p0); + FacetLocation f1 = distanceToB.nearestLocation(seg.p1); + if (f0.isSameSegment(f1)) + return true; + //-- check if endpoints of same segment + if (f1.getIndex() == f0.normalize(f0.getIndex() + 1) + && seg.p1.equals2D(f1.getEndPoint(0))) + { + return true; + } + if (f0.getIndex() == f0.normalize(f1.getIndex() + 1) + && seg.p0.equals2D(f0.getEndPoint(0))) + { + return true; + } + return false; + } + private static boolean isBeyondLimit(double maxDist, double maxDistanceLimit) { return maxDistanceLimit >= 0 && maxDist > maxDistanceLimit; } @@ -478,6 +506,10 @@ public TargetDistance(Geometry geom) { } } + public FacetLocation nearestLocation(Coordinate p) { + return distanceToFacets.nearestLocation(p); + } + public Coordinate nearestFacetPoint(Coordinate p) { return distanceToFacets.nearestPoint(p); } From 85e964ebb25401eedf5cc140517d385c367a4dd4 Mon Sep 17 00:00:00 2001 From: Martin Davis Date: Mon, 9 Feb 2026 15:59:14 -0800 Subject: [PATCH 20/43] Improve checking for collinear segments --- .../distance/DirectedHausdorffDistance.java | 18 +------ .../jts/operation/distance/FacetLocation.java | 52 ++++++++++++++++++- .../jts/operation/distance/FacetSequence.java | 12 ++--- 3 files changed, 58 insertions(+), 24 deletions(-) diff --git a/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java b/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java index 8758c826b8..81af952703 100644 --- a/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java +++ b/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java @@ -314,7 +314,7 @@ private Coordinate[] computeAtEdges(Geometry geomA, double tolerance, double max * - if segment distance bound is less than distance limit, no other segment can be farther * - if a point of segment is farther than limit, isFulyWithin must be false */ - if (segMaxDist.length() <= tolerance + if ((segMaxDist.length() <= tolerance) || isWithinLimit(segMaxDist.maxDistanceBound(), maxDistanceLimit) || isBeyondLimit(segMaxDist.maxDistance(), maxDistanceLimit) ) { @@ -328,7 +328,6 @@ private Coordinate[] computeAtEdges(Geometry geomA, double tolerance, double max continue; //System.out.println(segMaxDist); isSameSegment(segMaxDist); - } //System.out.println(segMaxDist); @@ -354,20 +353,7 @@ private Coordinate[] computeAtEdges(Geometry geomA, double tolerance, double max private boolean isSameSegment(DHDSegment seg) { FacetLocation f0 = distanceToB.nearestLocation(seg.p0); FacetLocation f1 = distanceToB.nearestLocation(seg.p1); - if (f0.isSameSegment(f1)) - return true; - //-- check if endpoints of same segment - if (f1.getIndex() == f0.normalize(f0.getIndex() + 1) - && seg.p1.equals2D(f1.getEndPoint(0))) - { - return true; - } - if (f0.getIndex() == f0.normalize(f1.getIndex() + 1) - && seg.p0.equals2D(f0.getEndPoint(0))) - { - return true; - } - return false; + return f0.isSameSegment(f1); } private static boolean isBeyondLimit(double maxDist, double maxDistanceLimit) { diff --git a/modules/core/src/main/java/org/locationtech/jts/operation/distance/FacetLocation.java b/modules/core/src/main/java/org/locationtech/jts/operation/distance/FacetLocation.java index 3b0d543b87..4331b8e6d8 100644 --- a/modules/core/src/main/java/org/locationtech/jts/operation/distance/FacetLocation.java +++ b/modules/core/src/main/java/org/locationtech/jts/operation/distance/FacetLocation.java @@ -1,9 +1,32 @@ +/* + * Copyright (c) 2026 Martin Davis + * + * All rights reserved. This program and the accompanying materials + * are made available under the terms of the Eclipse Public License 2.0 + * and Eclipse Distribution License v. 1.0 which accompanies this distribution. + * The Eclipse Public License is available at http://www.eclipse.org/legal/epl-v20.html + * and the Eclipse Distribution License is available at + * + * http://www.eclipse.org/org/documents/edl-v10.php. + */ package org.locationtech.jts.operation.distance; import org.locationtech.jts.geom.Coordinate; import org.locationtech.jts.geom.CoordinateSequence; import org.locationtech.jts.geom.CoordinateSequences; +/** + * + * Location indexes are always the index of a sequence segment. + * This means they are always less than the number of vertices + * in the sequence. The endpoint in a sequence + * has the index of the final segment in the sequence. + * if the sequence is a ring, the imdex of the final endpoint is + * normalized to 0. + * + * @author mdavis + * + */ public class FacetLocation { private CoordinateSequence seq; private int index; @@ -11,12 +34,37 @@ public class FacetLocation { public FacetLocation(CoordinateSequence seq, int index, Coordinate pt) { this.seq = seq; - this.index = index; this.pt = pt; + this.index = index; + if (index >= seq.size()) + this.index = seq.size() - 1; } public boolean isSameSegment(FacetLocation f) { - return seq == f.seq && index == f.index; + if (seq != f.seq) + return false; + if (index == f.index) + return true; + //-- check for end pt same as start point of next segment + if (isNext(index, f.index)) { + Coordinate endPt = seq.getCoordinate(index + 1); + return f.pt.equals2D(endPt); + } + if (isNext(f.index, index)) { + Coordinate endPt = f.seq.getCoordinate(index + 1); + return pt.equals2D(endPt); + } + return false; + } + + private boolean isNext(int index, int index1) { + if (index1 == index + 1) + return true; + if (index1 == 0 && CoordinateSequences.isRing(seq) + && index1 == seq.size() - 1) { + return true; + } + return false; } public int getIndex() { diff --git a/modules/core/src/main/java/org/locationtech/jts/operation/distance/FacetSequence.java b/modules/core/src/main/java/org/locationtech/jts/operation/distance/FacetSequence.java index ce07747f55..7b1171c59c 100644 --- a/modules/core/src/main/java/org/locationtech/jts/operation/distance/FacetSequence.java +++ b/modules/core/src/main/java/org/locationtech/jts/operation/distance/FacetSequence.java @@ -154,12 +154,9 @@ else if (isPointOther) { public FacetLocation nearestLocation(Coordinate p) { - boolean isPoint = isPoint(); - - if (isPoint) { + if (isPoint()) { return new FacetLocation(pts, 0, pts.getCoordinate(0)); } - return nearestLocationOnLine(p); } @@ -177,9 +174,12 @@ private FacetLocation nearestLocationOnLine(Coordinate pt) minDistance = dist; nearestPt = nearestPoint(pt, q0, q1); index = i; - //-- segments are half-open, so final endpoint belongs to next segment + //-- segments are half-open, so 2nd endpoint belongs to next segment + //-- except for last segment on non-closed sequence if (dist == 0.0 && pt.equals2D(q1)) { - index++; + if (index < pts.size() - 1) { + index++; + } //-- normalize index for a ring index = normalize(pts, index); } From c670eccd1ebc52367ed04993b31209b5a82f4037 Mon Sep 17 00:00:00 2001 From: Martin Davis Date: Mon, 9 Feb 2026 22:07:38 -0800 Subject: [PATCH 21/43] Improve segment search termination criteria --- .../distance/DirectedHausdorffDistance.java | 59 +++++++++++++------ .../DirectedHausdorffDistanceTest.java | 4 +- 2 files changed, 43 insertions(+), 20 deletions(-) diff --git a/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java b/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java index 81af952703..31b24593a3 100644 --- a/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java +++ b/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java @@ -299,45 +299,66 @@ private Coordinate[] computeAtPoints(Geometry geomA, double maxDistanceLimit) { private Coordinate[] computeAtEdges(Geometry geomA, double tolerance, double maxDistanceLimit) { PriorityQueue segQueue = createSegQueue(geomA); - boolean maxDistSegFound = false; DHDSegment segMaxDist = null; long iter = 0; while (! segQueue.isEmpty()) { iter++; - // get the segment with greatest distance - segMaxDist = segQueue.remove(); - + // get the segment with greatest distance bound + DHDSegment segMaxBound = segQueue.remove(); +/* + double maxDistBound = segMaxBound.maxDistanceBound(); + double maxDist = segMaxBound.maxDistance(); + System.out.format("%s len: %f bound: %f maxDist: %f\n", + segMaxBound, segMaxBound.length(), maxDistBound, maxDist); + if (maxDist > 1) { + System.out.println("FOUND"); + } +*/ + /** + * Save if segment is farther than most distant so far + */ + if (segMaxDist == null + || segMaxBound.maxDistance() > segMaxDist.maxDistance()) { + segMaxDist = segMaxBound; + } /** - * Exit if segment length is within tolerance, so distance is accurate enough. - * * If maxDistanceLimit is specified, short-circuit if: * - if segment distance bound is less than distance limit, no other segment can be farther * - if a point of segment is farther than limit, isFulyWithin must be false */ - if ((segMaxDist.length() <= tolerance) - || isWithinLimit(segMaxDist.maxDistanceBound(), maxDistanceLimit) - || isBeyondLimit(segMaxDist.maxDistance(), maxDistanceLimit) + if (isWithinLimit(segMaxBound.maxDistanceBound(), maxDistanceLimit) + || isBeyondLimit(segMaxBound.maxDistance(), maxDistanceLimit) ) { - maxDistSegFound = true; break; } - //-- check for equal or coincident segments - if (segMaxDist.maxDistance() == 0.0) { - if (isSameSegment(segMaxDist)) + /** + * Check for equal or coincident segments. + * If so, don't bisect the segment further. + * This improves performance when the inputs have identical segments. + */ + if (segMaxBound.maxDistance() == 0.0) { + if (isSameSegment(segMaxBound)) continue; //System.out.println(segMaxDist); - isSameSegment(segMaxDist); + //isSameSegment(segMaxBound); } //System.out.println(segMaxDist); - //-- not within tolerance, so bisect segment and keep searching - DHDSegment[] bisects = segMaxDist.bisect(distanceToB); - addNonInterior(bisects[0], segQueue); - addNonInterior(bisects[1], segQueue); + /** + * If segment is longer than tolerance + * and it might provide a better max distance point, + * bisect and keep searching + */ + if ((segMaxBound.length() > tolerance) + && segMaxBound.maxDistanceBound() > segMaxDist.maxDistance()) { + DHDSegment[] bisects = segMaxBound.bisect(distanceToB); + addNonInterior(bisects[0], segQueue); + addNonInterior(bisects[1], segQueue); + } } - if (maxDistSegFound) + if (segMaxDist != null) return segMaxDist.getMaxDistPts(); /** diff --git a/modules/core/src/test/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistanceTest.java b/modules/core/src/test/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistanceTest.java index 58e660111d..7e10145f2c 100644 --- a/modules/core/src/test/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistanceTest.java +++ b/modules/core/src/test/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistanceTest.java @@ -132,7 +132,9 @@ public void testMultiPolygons() String a = "MULTIPOLYGON (((1 1, 1 10, 5 1, 1 1)), ((4 17, 9 15, 9 6, 4 17)))"; String b = "MULTIPOLYGON (((1 12, 4 13, 8 10, 1 12)), ((3 8, 7 7, 6 2, 3 8)))"; checkDistance(a, b, 0.01, "LINESTRING (1 1, 5.4 3.2)"); - checkDistance(b, a, 0.001, "LINESTRING (2.669921875 12.556640625, 5.446115154109589 13.818546660958905)"); + checkDistanceStartPtLen(b, a, 0.001, + "LINESTRING (2.669921875 12.556640625, 5.446115154109589 13.818546660958905)", + 0.01); } // Tests that target area interiors have distance = 0 From ed9c34d817c71b5ec1503a8fc9249b86beeeba0d Mon Sep 17 00:00:00 2001 From: Martin Davis Date: Mon, 9 Feb 2026 22:19:16 -0800 Subject: [PATCH 22/43] Optimze adding segments to queue --- .../distance/DirectedHausdorffDistance.java | 26 ++++++++++++------- 1 file changed, 16 insertions(+), 10 deletions(-) diff --git a/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java b/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java index 31b24593a3..92d95c6bc4 100644 --- a/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java +++ b/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java @@ -80,11 +80,6 @@ * This algorithm is easier to use, more accurate, * and much faster than {@link DiscreteHausdorffDistance}. * - *

KNOWN ISSUES

- *
    - *
  • if the two geometries are identical or nearly so, - * performance may be slow. - *
* @author Martin Davis * */ @@ -481,6 +476,7 @@ public void filter(Geometry geom) { * @param priq */ private void addSegments(Coordinate[] pts, PriorityQueue priq) { + DHDSegment segMaxDist = null; DHDSegment prevSeg = null; for (int i = 0; i < pts.length - 1; i++) { DHDSegment seg; @@ -488,14 +484,24 @@ private void addSegments(Coordinate[] pts, PriorityQueue priq) { seg = DHDSegment.create(pts[i], pts[i + 1], distanceToB); } else { - //-- optimize by avoiding recomputing pt distance + //-- avoiding recomputing prev pt distance seg = DHDSegment.create(prevSeg, pts[i + 1], distanceToB); } prevSeg = seg; - /** - * Don't add interior segments, since their distance must be zero. - */ - addNonInterior(seg, priq); + + //-- don't bother adding segment if it can't be further away then current max + if (segMaxDist == null + || seg.maxDistanceBound() > segMaxDist.maxDistance()) { + /** + * Don't add interior segments, since their distance must be zero. + */ + addNonInterior(seg, priq); + } + + if (segMaxDist == null + || seg.maxDistance() > segMaxDist.maxDistance()) { + segMaxDist = seg;; + } //System.out.println(seg.distance()); } } From a02d44eecf8dd0f48c30d60fab5641e1dcc11f2f Mon Sep 17 00:00:00 2001 From: Martin Davis Date: Tue, 10 Feb 2026 08:49:42 -0800 Subject: [PATCH 23/43] cache maxDistance --- .../distance/DirectedHausdorffDistance.java | 36 +++++++++---------- 1 file changed, 17 insertions(+), 19 deletions(-) diff --git a/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java b/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java index 92d95c6bc4..f278e49fb2 100644 --- a/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java +++ b/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java @@ -328,15 +328,15 @@ private Coordinate[] computeAtEdges(Geometry geomA, double tolerance, double max } /** - * Check for equal or coincident segments. + * Check for equal or collinear segments. * If so, don't bisect the segment further. - * This improves performance when the inputs have identical segments. + * This greatly improves performance when the inputs + * have identical or collinear segments + * (in particular, the case when the inputs are identical). */ if (segMaxBound.maxDistance() == 0.0) { - if (isSameSegment(segMaxBound)) + if (isSameOrCollinear(segMaxBound)) continue; - //System.out.println(segMaxDist); - //isSameSegment(segMaxBound); } //System.out.println(segMaxDist); @@ -366,7 +366,7 @@ private Coordinate[] computeAtEdges(Geometry geomA, double tolerance, double max return pair(maxPt, maxPt); } - private boolean isSameSegment(DHDSegment seg) { + private boolean isSameOrCollinear(DHDSegment seg) { FacetLocation f0 = distanceToB.nearestLocation(seg.p0); FacetLocation f1 = distanceToB.nearestLocation(seg.p1); return f0.isSameSegment(f1); @@ -471,6 +471,7 @@ public void filter(Geometry geom) { } /** + * Add segments to queue * * @param pts * @param priq @@ -489,7 +490,7 @@ private void addSegments(Coordinate[] pts, PriorityQueue priq) { } prevSeg = seg; - //-- don't bother adding segment if it can't be further away then current max + //-- don't add segment if it can't be further away then current max if (segMaxDist == null || seg.maxDistanceBound() > segMaxDist.maxDistance()) { /** @@ -575,7 +576,8 @@ public static DHDSegment create(DHDSegment prevSeg, Coordinate p1, TargetDistanc private Coordinate nearPt0; private Coordinate p1; private Coordinate nearPt1; - double maxDistanceBound = Double.NEGATIVE_INFINITY; + private double maxDistanceBound = Double.NEGATIVE_INFINITY; + private double maxDistance; private DHDSegment(Coordinate p0, Coordinate p1) { this.p0 = p0; @@ -587,19 +589,19 @@ private DHDSegment(Coordinate p0, Coordinate nearPt0, Coordinate p1, Coordinate this.nearPt0 = nearPt0; this.p1 = p1; this.nearPt1 = nearPt1; - computeMaxDistanceBound(); + computeMaxDistances(); } private void init(TargetDistance dist) { nearPt0 = dist.nearestPoint(p0); nearPt1 = dist.nearestPoint(p1); - computeMaxDistanceBound(); + computeMaxDistances(); } private void init(Coordinate nearest0, TargetDistance dist) { nearPt0 = nearest0; nearPt1 = dist.nearestPoint(p1); - computeMaxDistanceBound(); + computeMaxDistances(); } public Coordinate getEndpoint(int index) { @@ -611,9 +613,7 @@ public double length() { } public double maxDistance() { - double dist0 = p0.distance(nearPt0); - double dist1 = p1.distance(nearPt1); - return Math.max(dist0,dist1); + return maxDistance; } public double maxDistanceBound() { @@ -632,17 +632,15 @@ public Coordinate[] getMaxDistPts() { /** * Computes a least upper bound for the maximum distance to a segment. */ - private void computeMaxDistanceBound() { - //System.out.println(distance()); - + private void computeMaxDistances() { /** * Least upper bound is the max distance to the endpoints, * plus half segment length. */ double dist0 = p0.distance(nearPt0); double dist1 = p1.distance(nearPt1); - double maxDist = Math.max(dist0, dist1); - maxDistanceBound = maxDist + length() / 2; + maxDistance = Math.max(dist0, dist1); + maxDistanceBound = maxDistance + length() / 2; } public DHDSegment[] bisect(TargetDistance dist) { From bdfcef182583a6660a23b6ed7a1cce0445c5d5ee Mon Sep 17 00:00:00 2001 From: Martin Davis Date: Tue, 10 Feb 2026 08:53:54 -0800 Subject: [PATCH 24/43] Rename method --- .../distance/DirectedHausdorffDistance.java | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java b/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java index f278e49fb2..da95328dcc 100644 --- a/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java +++ b/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java @@ -69,7 +69,7 @@ *

* The class can be used in prepared mode. * Creating an instance on a target geometry caches indexes for that geometry. - * Then {@link #computeDistancePoints(Geometry, double) + * Then {@link #maximumDistancePoints(Geometry, double) * or {@link #isFullyWithinDistance(Geometry, double, double)} * can be called efficiently for multiple query geometries. *

@@ -96,7 +96,7 @@ public class DirectedHausdorffDistance { public static double distance(Geometry a, Geometry b, double tolerance) { DirectedHausdorffDistance hd = new DirectedHausdorffDistance(b); - return distance(hd.computeDistancePoints(a, tolerance)); + return distance(hd.maximumDistancePoints(a, tolerance)); } /** @@ -111,7 +111,7 @@ public static double distance(Geometry a, Geometry b, double tolerance) public static LineString distanceLine(Geometry a, Geometry b, double tolerance) { DirectedHausdorffDistance hd = new DirectedHausdorffDistance(b); - return a.getFactory().createLineString(hd.computeDistancePoints(a, tolerance)); + return a.getFactory().createLineString(hd.maximumDistancePoints(a, tolerance)); } /** @@ -144,9 +144,9 @@ public static boolean isFullyWithinDistance(Geometry a, Geometry b, double maxDi public static LineString hausdorffDistanceLine(Geometry a, Geometry b, double tolerance) { DirectedHausdorffDistance hdAB = new DirectedHausdorffDistance(b); - Coordinate[] ptsAB = hdAB.computeDistancePoints(a, tolerance); + Coordinate[] ptsAB = hdAB.maximumDistancePoints(a, tolerance); DirectedHausdorffDistance hdBA = new DirectedHausdorffDistance(a); - Coordinate[] ptsBA = hdBA.computeDistancePoints(b, tolerance); + Coordinate[] ptsBA = hdBA.maximumDistancePoints(b, tolerance); //-- return points in A-B order Coordinate[] pts; @@ -227,13 +227,13 @@ private static boolean isBeyond(Envelope envA, Envelope envB, double maxDistance /** * Computes a pair of points which attain the directed Hausdorff distance - * of a query geometry A from the target. + * of a query geometry A from the target B. * * @param geomA the query geometry * @param tolerance the approximation distance tolerance * @return a pair of points [ptA, ptB] demonstrating the distance */ - public Coordinate[] computeDistancePoints(Geometry geomA, double tolerance) { + public Coordinate[] maximumDistancePoints(Geometry geomA, double tolerance) { return computeDistancePoints(geomA, tolerance, -1.0); } From 99c237ca3d292218ad34411c9ef4debe18e777e4 Mon Sep 17 00:00:00 2001 From: Martin Davis Date: Tue, 10 Feb 2026 15:08:11 -0800 Subject: [PATCH 25/43] Add topo equal lines test --- .../algorithm/distance/DirectedHausdorffDistanceTest.java | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/modules/core/src/test/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistanceTest.java b/modules/core/src/test/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistanceTest.java index 7e10145f2c..ab6c1880df 100644 --- a/modules/core/src/test/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistanceTest.java +++ b/modules/core/src/test/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistanceTest.java @@ -77,6 +77,14 @@ public void testLinePoints() "LINESTRING (0 0, 0 2)"); } + public void testLinesTopoEqual() + { + checkDistance("MULTILINESTRING ((10 10, 10 90, 40 30), (40 30, 60 80, 90 30, 40 10))", + "LINESTRING (10 10, 10 90, 40 30, 60 80, 90 30, 40 10)", + 0.01, + 0.0); + } + public void testLinesPolygon() { checkHD("MULTILINESTRING ((1 1, 2 7), (7 1, 9 9))", From ab16df31e77d84c45f76af71fd1b6676b3cda30a Mon Sep 17 00:00:00 2001 From: Martin Davis Date: Thu, 12 Feb 2026 19:29:03 -0800 Subject: [PATCH 26/43] Rename CoordinateSequenceLocation --- .../distance/DirectedHausdorffDistance.java | 8 ++++---- ...ation.java => CoordinateSequenceLocation.java} | 15 ++++++++++----- .../jts/operation/distance/FacetSequence.java | 8 ++++---- .../operation/distance/IndexedFacetDistance.java | 4 ++-- 4 files changed, 20 insertions(+), 15 deletions(-) rename modules/core/src/main/java/org/locationtech/jts/operation/distance/{FacetLocation.java => CoordinateSequenceLocation.java} (83%) diff --git a/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java b/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java index da95328dcc..151faaefbb 100644 --- a/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java +++ b/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java @@ -26,7 +26,7 @@ import org.locationtech.jts.geom.Location; import org.locationtech.jts.geom.Point; import org.locationtech.jts.io.WKTWriter; -import org.locationtech.jts.operation.distance.FacetLocation; +import org.locationtech.jts.operation.distance.CoordinateSequenceLocation; import org.locationtech.jts.operation.distance.IndexedFacetDistance; /** @@ -367,8 +367,8 @@ private Coordinate[] computeAtEdges(Geometry geomA, double tolerance, double max } private boolean isSameOrCollinear(DHDSegment seg) { - FacetLocation f0 = distanceToB.nearestLocation(seg.p0); - FacetLocation f1 = distanceToB.nearestLocation(seg.p1); + CoordinateSequenceLocation f0 = distanceToB.nearestLocation(seg.p0); + CoordinateSequenceLocation f1 = distanceToB.nearestLocation(seg.p1); return f0.isSameSegment(f1); } @@ -520,7 +520,7 @@ public TargetDistance(Geometry geom) { } } - public FacetLocation nearestLocation(Coordinate p) { + public CoordinateSequenceLocation nearestLocation(Coordinate p) { return distanceToFacets.nearestLocation(p); } diff --git a/modules/core/src/main/java/org/locationtech/jts/operation/distance/FacetLocation.java b/modules/core/src/main/java/org/locationtech/jts/operation/distance/CoordinateSequenceLocation.java similarity index 83% rename from modules/core/src/main/java/org/locationtech/jts/operation/distance/FacetLocation.java rename to modules/core/src/main/java/org/locationtech/jts/operation/distance/CoordinateSequenceLocation.java index 4331b8e6d8..d0d03c5a1b 100644 --- a/modules/core/src/main/java/org/locationtech/jts/operation/distance/FacetLocation.java +++ b/modules/core/src/main/java/org/locationtech/jts/operation/distance/CoordinateSequenceLocation.java @@ -16,23 +16,24 @@ import org.locationtech.jts.geom.CoordinateSequences; /** + * A location on a {@link FacetSequence}. * * Location indexes are always the index of a sequence segment. - * This means they are always less than the number of vertices + * Thus they are always less than the number of vertices * in the sequence. The endpoint in a sequence * has the index of the final segment in the sequence. - * if the sequence is a ring, the imdex of the final endpoint is + * if the sequence is a ring, the index of the final endpoint is * normalized to 0. * * @author mdavis * */ -public class FacetLocation { +public class CoordinateSequenceLocation { private CoordinateSequence seq; private int index; private Coordinate pt; - public FacetLocation(CoordinateSequence seq, int index, Coordinate pt) { + public CoordinateSequenceLocation(CoordinateSequence seq, int index, Coordinate pt) { this.seq = seq; this.pt = pt; this.index = index; @@ -40,7 +41,11 @@ public FacetLocation(CoordinateSequence seq, int index, Coordinate pt) { this.index = seq.size() - 1; } - public boolean isSameSegment(FacetLocation f) { + public Coordinate getCoordinate() { + return pt; + } + + public boolean isSameSegment(CoordinateSequenceLocation f) { if (seq != f.seq) return false; if (index == f.index) diff --git a/modules/core/src/main/java/org/locationtech/jts/operation/distance/FacetSequence.java b/modules/core/src/main/java/org/locationtech/jts/operation/distance/FacetSequence.java index 7b1171c59c..efe70c3b0d 100644 --- a/modules/core/src/main/java/org/locationtech/jts/operation/distance/FacetSequence.java +++ b/modules/core/src/main/java/org/locationtech/jts/operation/distance/FacetSequence.java @@ -152,15 +152,15 @@ else if (isPointOther) { return locs; } - public FacetLocation nearestLocation(Coordinate p) + public CoordinateSequenceLocation nearestLocation(Coordinate p) { if (isPoint()) { - return new FacetLocation(pts, 0, pts.getCoordinate(0)); + return new CoordinateSequenceLocation(pts, 0, pts.getCoordinate(0)); } return nearestLocationOnLine(p); } - private FacetLocation nearestLocationOnLine(Coordinate pt) + private CoordinateSequenceLocation nearestLocationOnLine(Coordinate pt) { double minDistance = Double.MAX_VALUE; int index = -1; @@ -187,7 +187,7 @@ private FacetLocation nearestLocationOnLine(Coordinate pt) break; } } - return new FacetLocation(pts, index, nearestPt); + return new CoordinateSequenceLocation(pts, index, nearestPt); } private static int normalize(CoordinateSequence pts, int index) { diff --git a/modules/core/src/main/java/org/locationtech/jts/operation/distance/IndexedFacetDistance.java b/modules/core/src/main/java/org/locationtech/jts/operation/distance/IndexedFacetDistance.java index 3615ca6ae3..bd8ff34ea4 100644 --- a/modules/core/src/main/java/org/locationtech/jts/operation/distance/IndexedFacetDistance.java +++ b/modules/core/src/main/java/org/locationtech/jts/operation/distance/IndexedFacetDistance.java @@ -139,7 +139,7 @@ public Coordinate[] nearestPoints(Geometry g) return fs1.nearestLocations(fs2); } - public FacetLocation nearestLocation(Coordinate p) { + public CoordinateSequenceLocation nearestLocation(Coordinate p) { CoordinateSequence seq = new CoordinateArraySequence(new Coordinate[] { p }); FacetSequence fs = new FacetSequence(seq, 0); Object nearest = cachedTree.nearestNeighbour(fs.getEnvelope(), fs, FACET_SEQ_DIST); @@ -161,7 +161,7 @@ public Coordinate nearestPoint(Coordinate p) FacetSequence fs = new FacetSequence(seq, 0); Object nearest = cachedTree.nearestNeighbour(fs.getEnvelope(), fs, FACET_SEQ_DIST); FacetSequence fsN = (FacetSequence) nearest; - return fsN.nearestLocations(fs)[0]; + return fsN.nearestLocation(p).getCoordinate(); } /** From 19628ca1a4b578ddac0269df45f6820ed66cef0f Mon Sep 17 00:00:00 2001 From: Martin Davis Date: Fri, 20 Feb 2026 16:25:01 -0800 Subject: [PATCH 27/43] Rename distance functions --- .../jtstest/function/DistanceFunctions.java | 94 ++++++++----------- 1 file changed, 39 insertions(+), 55 deletions(-) diff --git a/modules/app/src/main/java/org/locationtech/jtstest/function/DistanceFunctions.java b/modules/app/src/main/java/org/locationtech/jtstest/function/DistanceFunctions.java index 2926529470..65c74f0a5e 100644 --- a/modules/app/src/main/java/org/locationtech/jtstest/function/DistanceFunctions.java +++ b/modules/app/src/main/java/org/locationtech/jtstest/function/DistanceFunctions.java @@ -46,53 +46,60 @@ public static Geometry frechetDistanceLine(Geometry a, Geometry b) return a.getFactory().createLineString(dist.getCoordinates()); } - public static double hausdorffDistance(Geometry a, Geometry b) + @Metadata(description="Oriented discrete Hausdorff distance line from A to B") + public static Geometry orientedDiscreteHausdorffLine(Geometry a, Geometry b) { - return DiscreteHausdorffDistance.distance(a, b); - } - - @Metadata(description="Hausdorff distance between A and B") - public static Geometry hausdorffDistanceLine(Geometry a, Geometry b) - { - return DiscreteHausdorffDistance.distanceLine(a, b); + return DiscreteHausdorffDistance.orientedDistanceLine(a, b); } - @Metadata(description="Hausdorff distance between A and B, densified") - public static Geometry hausdorffDistanceLineDensify(Geometry a, Geometry b, - @Metadata(title="Densify fraction") - double frac) + @Metadata(description="Oriented discrete Hausdorff distance from A to B") + public static double orientedDiscreteHausdorffDistance(Geometry a, Geometry b) { - return DiscreteHausdorffDistance.distanceLine(a, b, frac); + return DiscreteHausdorffDistance.orientedDistance(a, b); } - - @Metadata(description="Oriented Hausdorff distance from A to B") - public static Geometry orientedHausdorffDistanceLine(Geometry a, Geometry b) + + @Metadata(description="Oriented discrete Hausdorff distance line from A to B, densified") + public static Geometry orientedDiscreteHausdorffLineDensify(Geometry a, Geometry b, + @Metadata(title="Densify fraction") + double frac) { - return DiscreteHausdorffDistance.orientedDistanceLine(a, b); + return DiscreteHausdorffDistance.orientedDistanceLine(a, b, frac); } - @Metadata(description="Oriented Hausdorff distance from A to B") + @Metadata(description="Clipped Oriented Hausdorff distance from A to B") public static Geometry clippedOrientedHausdorffDistanceLine(Geometry a, Geometry b) { //TODO: would this be more efficient done as part of DiscreteHausdorffDistance? Geometry clippedLine = LinearReferencingFunctions.project(a, b); return DiscreteHausdorffDistance.orientedDistanceLine(clippedLine, b); } - - @Metadata(description="Oriented Hausdorff distance from A to B") - public static double orientedHausdorffDistance(Geometry a, Geometry b) - { - return DiscreteHausdorffDistance.orientedDistance(a, b); - } - - @Metadata(description="Oriented Hausdorff distance from A to B, densified") - public static Geometry orientedHausdorffDistanceLineDensify(Geometry a, Geometry b, - @Metadata(title="Densify fraction") - double frac) + + @Metadata(description="Directed Hausdorff distance from A to B, up to tolerance") + public static double directedHausdorffDistance(Geometry a, Geometry b, + @Metadata(title="Distance tolerance") + double distTol) { - return DiscreteHausdorffDistance.orientedDistanceLine(a, b, frac); + return DirectedHausdorffDistance.distance(a, b, distTol); } - + + @Metadata(description="Directed Hausdorff distance from A to B, up to tolerance") + public static Geometry directedHausdorffLine(Geometry a, Geometry b, + @Metadata(title="Distance tolerance") + double distTol) + { + return DirectedHausdorffDistance.distanceLine(a, b, distTol); + } + + @Metadata(description="Hausdorff distance between A and B, up to tolerance") + public static Geometry hausdorffLine(Geometry a, Geometry b, + @Metadata(title="Distance tolerance") + double distTol) + { + return DirectedHausdorffDistance.hausdorffDistanceLine(a, b, distTol); + } + + //-------------------------------------------- + public static double distanceIndexed(Geometry a, Geometry b) { return IndexedFacetDistance.distance(a, b); } @@ -118,28 +125,5 @@ public static Geometry nearestPointsIndexedEachB(Geometry a, Geometry b) { return a.getFactory().createMultiLineString(lines); } - - @Metadata(description="Directed Hausdorff distance from A to B, up to tolerance") - public static double dhd(Geometry a, Geometry b, - @Metadata(title="Distance tolerance") - double distTol) - { - return DirectedHausdorffDistance.distance(a, b, distTol); - } - - @Metadata(description="Directed Hausdorff distance from A to B, up to tolerance") - public static Geometry dhdLine(Geometry a, Geometry b, - @Metadata(title="Distance tolerance") - double distTol) - { - return DirectedHausdorffDistance.distanceLine(a, b, distTol); - } - - @Metadata(description="Hausdorff distance between A and B, up to tolerance") - public static Geometry hdLine(Geometry a, Geometry b, - @Metadata(title="Distance tolerance") - double distTol) - { - return DirectedHausdorffDistance.hausdorffDistanceLine(a, b, distTol); - } + } From 3705cc1c430ef98d5b29f77a5b77ba0c192a5862 Mon Sep 17 00:00:00 2001 From: Martin Davis Date: Fri, 20 Feb 2026 16:35:52 -0800 Subject: [PATCH 28/43] Method renames --- .../distance/DirectedHausdorffDistance.java | 82 ++++++++++--------- 1 file changed, 43 insertions(+), 39 deletions(-) diff --git a/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java b/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java index 151faaefbb..973de87f53 100644 --- a/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java +++ b/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java @@ -100,7 +100,7 @@ public static double distance(Geometry a, Geometry b, double tolerance) } /** - * Computes a pair of points which attain the directed Hausdorff distance + * Computes a line containing a pair of points which attain the directed Hausdorff distance * of a query geometry A from a target one B. * * @param a the query geometry @@ -113,27 +113,9 @@ public static LineString distanceLine(Geometry a, Geometry b, double tolerance) DirectedHausdorffDistance hd = new DirectedHausdorffDistance(b); return a.getFactory().createLineString(hd.maximumDistancePoints(a, tolerance)); } - - /** - * Computes whether a query geometry lies fully within a give distance of a target geometry. - * Equivalently, detects whether any point of the query geometry is farther - * from the target than the specified distance. - * This is the case if DHD(A, B) > maxDistance. - * - * @param a the query geometry - * @param b the target geometry - * @param maxDistance the distance limit - * @param tolerance the approximation distance tolerance - * @return true if the query geometry lies fully within the distance of the target - */ - public static boolean isFullyWithinDistance(Geometry a, Geometry b, double maxDistance, double tolerance) - { - DirectedHausdorffDistance hd = new DirectedHausdorffDistance(b); - return hd.isFullyWithinDistance(a, maxDistance, tolerance); - } /** - * Computes a pair of points which attain the Hausdorff distance + * Computes a pair of points which attain the symmetric Hausdorff distance * between two geometries. * * @param a a geometry @@ -159,6 +141,24 @@ public static LineString hausdorffDistanceLine(Geometry a, Geometry b, double to return a.getFactory().createLineString(pts); } + /** + * Computes whether a query geometry lies fully within a give distance of a target geometry. + * Equivalently, detects whether any point of the query geometry is farther + * from the target than the specified distance. + * This is the case if DHD(A, B) > maxDistance. + * + * @param a the query geometry + * @param b the target geometry + * @param maxDistance the distance limit + * @param tolerance the approximation distance tolerance + * @return true if the query geometry lies fully within the distance of the target + */ + public static boolean isFullyWithinDistance(Geometry a, Geometry b, double maxDistance, double tolerance) + { + DirectedHausdorffDistance hd = new DirectedHausdorffDistance(b); + return hd.isFullyWithinDistance(a, maxDistance, tolerance); + } + private static double distance(Coordinate[] pts) { return pts[0].distance(pts[1]); } @@ -239,10 +239,10 @@ public Coordinate[] maximumDistancePoints(Geometry geomA, double tolerance) { private Coordinate[] computeDistancePoints(Geometry geomA, double tolerance, double maxDistanceLimit) { if (geomA.getDimension() == Dimension.P) { - return computeAtPoints(geomA, maxDistanceLimit); + return computeForPoints(geomA, maxDistanceLimit); } //TODO: handle mixed geoms with points - Coordinate[] maxDistPtsEdge = computeAtEdges(geomA, tolerance, maxDistanceLimit); + Coordinate[] maxDistPtsEdge = computeForEdges(geomA, tolerance, maxDistanceLimit); if (isBeyondLimit(distance(maxDistPtsEdge), maxDistanceLimit)) { return maxDistPtsEdge; @@ -252,7 +252,7 @@ private Coordinate[] computeDistancePoints(Geometry geomA, double tolerance, dou * Polygonal query geometry may have an interior point as the farthest point. */ if (geomA.getDimension() == Dimension.A) { - Coordinate[] maxDistPtsInterior = computeAtAreaInterior(geomA, tolerance); + Coordinate[] maxDistPtsInterior = computeForAreaInterior(geomA, tolerance); if (maxDistPtsInterior != null && distance(maxDistPtsInterior) > distance(maxDistPtsEdge)) { return maxDistPtsInterior; @@ -261,7 +261,7 @@ && distance(maxDistPtsInterior) > distance(maxDistPtsEdge)) { return maxDistPtsEdge; } - private Coordinate[] computeAtPoints(Geometry geomA, double maxDistanceLimit) { + private Coordinate[] computeForPoints(Geometry geomA, double maxDistanceLimit) { double maxDist = -1.0;; Coordinate[] maxDistPtsAB = null; Iterator geomi = new GeometryCollectionIterator(geomA); @@ -291,7 +291,7 @@ private Coordinate[] computeAtPoints(Geometry geomA, double maxDistanceLimit) { return maxDistPtsAB; } - private Coordinate[] computeAtEdges(Geometry geomA, double tolerance, double maxDistanceLimit) { + private Coordinate[] computeForEdges(Geometry geomA, double tolerance, double maxDistanceLimit) { PriorityQueue segQueue = createSegQueue(geomA); DHDSegment segMaxDist = null; @@ -313,7 +313,7 @@ private Coordinate[] computeAtEdges(Geometry geomA, double tolerance, double max * Save if segment is farther than most distant so far */ if (segMaxDist == null - || segMaxBound.maxDistance() > segMaxDist.maxDistance()) { + || segMaxBound.getMaxDistance() > segMaxDist.getMaxDistance()) { segMaxDist = segMaxBound; } /** @@ -321,8 +321,8 @@ private Coordinate[] computeAtEdges(Geometry geomA, double tolerance, double max * - if segment distance bound is less than distance limit, no other segment can be farther * - if a point of segment is farther than limit, isFulyWithin must be false */ - if (isWithinLimit(segMaxBound.maxDistanceBound(), maxDistanceLimit) - || isBeyondLimit(segMaxBound.maxDistance(), maxDistanceLimit) + if (isWithinLimit(segMaxBound.getMaxDistanceBound(), maxDistanceLimit) + || isBeyondLimit(segMaxBound.getMaxDistance(), maxDistanceLimit) ) { break; } @@ -334,7 +334,7 @@ private Coordinate[] computeAtEdges(Geometry geomA, double tolerance, double max * have identical or collinear segments * (in particular, the case when the inputs are identical). */ - if (segMaxBound.maxDistance() == 0.0) { + if (segMaxBound.getMaxDistance() == 0.0) { if (isSameOrCollinear(segMaxBound)) continue; } @@ -346,13 +346,17 @@ private Coordinate[] computeAtEdges(Geometry geomA, double tolerance, double max * and it might provide a better max distance point, * bisect and keep searching */ - if ((segMaxBound.length() > tolerance) - && segMaxBound.maxDistanceBound() > segMaxDist.maxDistance()) { + if ((segMaxBound.getLength() > tolerance) + && segMaxBound.getMaxDistanceBound() > segMaxDist.getMaxDistance()) { DHDSegment[] bisects = segMaxBound.bisect(distanceToB); addNonInterior(bisects[0], segQueue); addNonInterior(bisects[1], segQueue); } } + /** + * A segment at maximum distance was found. + * Return the farthest point pair + */ if (segMaxDist != null) return segMaxDist.getMaxDistPts(); @@ -395,7 +399,7 @@ private void addNonInterior(DHDSegment segment, PriorityQueue segQue * @return */ private boolean isInterior(DHDSegment segment) { - if (segment.maxDistance() > 0.0) { + if (segment.getMaxDistance() > 0.0) { return false; } return distanceToB.isInterior(segment.getEndpoint(0), segment.getEndpoint(1)); @@ -418,7 +422,7 @@ private boolean isInterior(DHDSegment segment) { * @return the maximum distance point pair at an interior point of A, * or null if it is known to not occur at an interior point */ - private Coordinate[] computeAtAreaInterior(Geometry geomA, double tolerance) { + private Coordinate[] computeForAreaInterior(Geometry geomA, double tolerance) { //TODO: extract polygonal geoms from A Geometry polygonalA = geomA; @@ -492,7 +496,7 @@ private void addSegments(Coordinate[] pts, PriorityQueue priq) { //-- don't add segment if it can't be further away then current max if (segMaxDist == null - || seg.maxDistanceBound() > segMaxDist.maxDistance()) { + || seg.getMaxDistanceBound() > segMaxDist.getMaxDistance()) { /** * Don't add interior segments, since their distance must be zero. */ @@ -500,7 +504,7 @@ private void addSegments(Coordinate[] pts, PriorityQueue priq) { } if (segMaxDist == null - || seg.maxDistance() > segMaxDist.maxDistance()) { + || seg.getMaxDistance() > segMaxDist.getMaxDistance()) { segMaxDist = seg;; } //System.out.println(seg.distance()); @@ -608,15 +612,15 @@ public Coordinate getEndpoint(int index) { return index == 0 ? p0 : p1; } - public double length() { + public double getLength() { return p0.distance(p1); } - public double maxDistance() { + public double getMaxDistance() { return maxDistance; } - public double maxDistanceBound() { + public double getMaxDistanceBound() { return maxDistanceBound; } @@ -640,7 +644,7 @@ private void computeMaxDistances() { double dist0 = p0.distance(nearPt0); double dist1 = p1.distance(nearPt1); maxDistance = Math.max(dist0, dist1); - maxDistanceBound = maxDistance + length() / 2; + maxDistanceBound = maxDistance + getLength() / 2; } public DHDSegment[] bisect(TargetDistance dist) { From 82f081a1158c5e1521f1e928c33438cf58883650 Mon Sep 17 00:00:00 2001 From: Martin Davis Date: Fri, 20 Feb 2026 16:39:49 -0800 Subject: [PATCH 29/43] Simplify code --- .../distance/DirectedHausdorffDistance.java | 15 +++++---------- 1 file changed, 5 insertions(+), 10 deletions(-) diff --git a/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java b/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java index 973de87f53..7b51139d46 100644 --- a/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java +++ b/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java @@ -117,6 +117,7 @@ public static LineString distanceLine(Geometry a, Geometry b, double tolerance) /** * Computes a pair of points which attain the symmetric Hausdorff distance * between two geometries. + * This the maximum of the two directed Hausdorff distances. * * @param a a geometry * @param b a geometry @@ -131,12 +132,10 @@ public static LineString hausdorffDistanceLine(Geometry a, Geometry b, double to Coordinate[] ptsBA = hdBA.maximumDistancePoints(b, tolerance); //-- return points in A-B order - Coordinate[] pts; - if (distance(ptsAB) > distance(ptsBA)) { - pts = ptsAB; - } - else { - pts = pairReverse(ptsBA); + Coordinate[] pts = ptsAB; + if (distance(ptsBA) > distance(ptsAB)) { + //-- reverse the BA points + pts = pair(ptsBA[1], ptsBA[0]); } return a.getFactory().createLineString(pts); } @@ -166,10 +165,6 @@ private static double distance(Coordinate[] pts) { private static Coordinate[] pair(Coordinate p0, Coordinate p1) { return new Coordinate[] { p0.copy(), p1.copy() }; } - - private static Coordinate[] pairReverse(Coordinate[] pts) { - return new Coordinate[] { pts[1].copy(), pts[0].copy() }; - } private Geometry geomB; private TargetDistance distanceToB; From 7685a90d5358f322dd469f0ca6eead07424750fa Mon Sep 17 00:00:00 2001 From: Martin Davis Date: Mon, 23 Feb 2026 11:16:38 -0800 Subject: [PATCH 30/43] Add auto-tolerance --- .../jtstest/function/DistanceFunctions.java | 10 +- .../jtstest/function/SelectionFunctions.java | 6 +- .../distance/DirectedHausdorffDistance.java | 128 ++++++++++++++++-- .../DirectedHausdorffDistanceTest.java | 2 +- 4 files changed, 128 insertions(+), 18 deletions(-) diff --git a/modules/app/src/main/java/org/locationtech/jtstest/function/DistanceFunctions.java b/modules/app/src/main/java/org/locationtech/jtstest/function/DistanceFunctions.java index 65c74f0a5e..71e6a177df 100644 --- a/modules/app/src/main/java/org/locationtech/jtstest/function/DistanceFunctions.java +++ b/modules/app/src/main/java/org/locationtech/jtstest/function/DistanceFunctions.java @@ -82,14 +82,20 @@ public static double directedHausdorffDistance(Geometry a, Geometry b, return DirectedHausdorffDistance.distance(a, b, distTol); } - @Metadata(description="Directed Hausdorff distance from A to B, up to tolerance") - public static Geometry directedHausdorffLine(Geometry a, Geometry b, + @Metadata(description="Directed Hausdorff distance line from A to B, up to tolerance") + public static Geometry directedHausdorffLineTol(Geometry a, Geometry b, @Metadata(title="Distance tolerance") double distTol) { return DirectedHausdorffDistance.distanceLine(a, b, distTol); } + @Metadata(description="Directed Hausdorff distance line from A to B") + public static Geometry directedHausdorffLine(Geometry a, Geometry b) + { + return DirectedHausdorffDistance.distanceLine(a, b); + } + @Metadata(description="Hausdorff distance between A and B, up to tolerance") public static Geometry hausdorffLine(Geometry a, Geometry b, @Metadata(title="Distance tolerance") diff --git a/modules/app/src/main/java/org/locationtech/jtstest/function/SelectionFunctions.java b/modules/app/src/main/java/org/locationtech/jtstest/function/SelectionFunctions.java index cc212f9ca9..799544c807 100644 --- a/modules/app/src/main/java/org/locationtech/jtstest/function/SelectionFunctions.java +++ b/modules/app/src/main/java/org/locationtech/jtstest/function/SelectionFunctions.java @@ -243,8 +243,7 @@ public static Geometry fullyWithinDistance(Geometry a, final Geometry mask, doub { return select(a, new GeometryPredicate() { public boolean isTrue(Geometry g) { - double tolerance = maximumDistance / 1e4; - return DirectedHausdorffDistance.isFullyWithinDistance(g, mask, maximumDistance, tolerance); + return DirectedHausdorffDistance.isFullyWithinDistance(g, mask, maximumDistance); } }); } @@ -254,8 +253,7 @@ public static Geometry fullyWithinDistancePrep(Geometry a, final Geometry mask, DirectedHausdorffDistance dhd = new DirectedHausdorffDistance(mask); return select(a, new GeometryPredicate() { public boolean isTrue(Geometry g) { - double tolerance = maximumDistance / 1e4; - return dhd.isFullyWithinDistance(g, maximumDistance, tolerance); + return dhd.isFullyWithinDistance(g, maximumDistance); } }); } diff --git a/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java b/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java index 7b51139d46..d8ea1344c4 100644 --- a/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java +++ b/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java @@ -86,11 +86,13 @@ public class DirectedHausdorffDistance { /** - * Computes the directed Hausdorff distance of a query geometry A from a target one B. + * Computes the directed Hausdorff distance + * of a query geometry A from a target one B, + * up to a given distance accuracy. * * @param a the query geometry * @param b the target geometry - * @param tolerance the approximation distance tolerance + * @param tolerance the accuracy distance tolerance * @return the directed Hausdorff distance */ public static double distance(Geometry a, Geometry b, double tolerance) @@ -100,12 +102,26 @@ public static double distance(Geometry a, Geometry b, double tolerance) } /** - * Computes a line containing a pair of points which attain the directed Hausdorff distance + * Computes the directed Hausdorff distance * of a query geometry A from a target one B. * * @param a the query geometry * @param b the target geometry - * @param tolerance the approximation distance tolerance + * @return the directed Hausdorff distance + */ + public static double distance(Geometry a, Geometry b) + { + DirectedHausdorffDistance hd = new DirectedHausdorffDistance(b); + return distance(hd.maximumDistancePoints(a)); + } + + /** + * Computes a line containing a pair of points which attain the directed Hausdorff distance + * of a query geometry A from a target one B, up to a given distance accuracy. + * + * @param a the query geometry + * @param b the target geometry + * @param tolerance the accuracy distance tolerance * @return a pair of points [ptA, ptB] demonstrating the distance */ public static LineString distanceLine(Geometry a, Geometry b, double tolerance) @@ -114,9 +130,24 @@ public static LineString distanceLine(Geometry a, Geometry b, double tolerance) return a.getFactory().createLineString(hd.maximumDistancePoints(a, tolerance)); } + /** + * Computes a line containing a pair of points which attain the directed Hausdorff distance + * of a query geometry A from a target one B. + * + * @param a the query geometry + * @param b the target geometry + * @param tolerance the accuracy distance tolerance + * @return a pair of points [ptA, ptB] demonstrating the distance + */ + public static LineString distanceLine(Geometry a, Geometry b) + { + DirectedHausdorffDistance hd = new DirectedHausdorffDistance(b); + return a.getFactory().createLineString(hd.maximumDistancePoints(a)); + } + /** * Computes a pair of points which attain the symmetric Hausdorff distance - * between two geometries. + * between two geometries, up to a given distance accuracy. * This the maximum of the two directed Hausdorff distances. * * @param a a geometry @@ -140,6 +171,31 @@ public static LineString hausdorffDistanceLine(Geometry a, Geometry b, double to return a.getFactory().createLineString(pts); } + /** + * Computes a pair of points which attain the symmetric Hausdorff distance + * between two geometries. + * This the maximum of the two directed Hausdorff distances. + * + * @param a a geometry + * @param b a geometry + * @return a pair of points [ptA, ptB] demonstrating the Hausdorff distance + */ + public static LineString hausdorffDistanceLine(Geometry a, Geometry b) + { + DirectedHausdorffDistance hdAB = new DirectedHausdorffDistance(b); + Coordinate[] ptsAB = hdAB.maximumDistancePoints(a); + DirectedHausdorffDistance hdBA = new DirectedHausdorffDistance(a); + Coordinate[] ptsBA = hdBA.maximumDistancePoints(b); + + //-- return points in A-B order + Coordinate[] pts = ptsAB; + if (distance(ptsBA) > distance(ptsAB)) { + //-- reverse the BA points + pts = pair(ptsBA[1], ptsBA[0]); + } + return a.getFactory().createLineString(pts); + } + /** * Computes whether a query geometry lies fully within a give distance of a target geometry. * Equivalently, detects whether any point of the query geometry is farther @@ -158,6 +214,12 @@ public static boolean isFullyWithinDistance(Geometry a, Geometry b, double maxDi return hd.isFullyWithinDistance(a, maxDistance, tolerance); } + public static boolean isFullyWithinDistance(Geometry a, Geometry b, double maxDistance) + { + DirectedHausdorffDistance hd = new DirectedHausdorffDistance(b); + return hd.isFullyWithinDistance(a, maxDistance); + } + private static double distance(Coordinate[] pts) { return pts[0].distance(pts[1]); } @@ -166,6 +228,20 @@ private static Coordinate[] pair(Coordinate p0, Coordinate p1) { return new Coordinate[] { p0.copy(), p1.copy() }; } + /** + * Heuristic factor to improve performance of area interior farthest point computation. + */ + private static final double AREA_INTERIOR_PERFORMANCE_FACTOR = 20; + + /** + * Heuristic automatic tolerance factor + */ + private static final double AUTO_TOLERANCE_FACTOR = 1.0e3; + + private static double computeTolerance(Geometry geom) { + return geom.getEnvelopeInternal().getDiameter() / AUTO_TOLERANCE_FACTOR; + } + private Geometry geomB; private TargetDistance distanceToB; @@ -180,14 +256,29 @@ public DirectedHausdorffDistance(Geometry b) { } /** - * Computes whether a query geometry lies fully within a give distance of the target geometry. + * Tests whether a query geometry lies fully within a give distance of the target geometry. * Equivalently, detects whether any point of the query geometry is farther * from the target than the specified distance. * This is the case if DHD(A, B) > maxDistance. * * @param a the query geometry * @param maxDistance the distance limit - * @param tolerance the approximation distance tolerance + * @return true if the query geometry lies fully within the distance of the target + */ + public boolean isFullyWithinDistance(Geometry a, double maxDistance) { + return isFullyWithinDistance(a, maxDistance, computeTolerance(a)); + } + + /** + * Tests whether a query geometry lies fully within a give distance of the target geometry, + * up to a given distance accuracy. + * Equivalently, detects whether any point of the query geometry is farther + * from the target than the specified distance. + * This is the case if DHD(A, B) > maxDistance. + * + * @param a the query geometry + * @param maxDistance the distance limit + * @param tolerance the accuracy distance tolerance * @return true if the query geometry lies fully within the distance of the target */ public boolean isFullyWithinDistance(Geometry a, double maxDistance, double tolerance) { @@ -219,7 +310,12 @@ private static boolean isBeyond(Envelope envA, Envelope envB, double maxDistance || envA.getMaxX() > envB.getMaxX() + maxDistance || envA.getMaxY() > envB.getMaxY() + maxDistance; } - + + public Coordinate[] maximumDistancePoints(Geometry geomA) { + double tolerance = computeTolerance(geomA); + return computeDistancePoints(geomA, tolerance, -1.0); + } + /** * Computes a pair of points which attain the directed Hausdorff distance * of a query geometry A from the target B. @@ -231,7 +327,7 @@ private static boolean isBeyond(Envelope envA, Envelope envB, double maxDistance public Coordinate[] maximumDistancePoints(Geometry geomA, double tolerance) { return computeDistancePoints(geomA, tolerance, -1.0); } - + private Coordinate[] computeDistancePoints(Geometry geomA, double tolerance, double maxDistanceLimit) { if (geomA.getDimension() == Dimension.P) { return computeForPoints(geomA, maxDistanceLimit); @@ -246,6 +342,7 @@ private Coordinate[] computeDistancePoints(Geometry geomA, double tolerance, dou /** * Polygonal query geometry may have an interior point as the farthest point. */ + //* if (geomA.getDimension() == Dimension.A) { Coordinate[] maxDistPtsInterior = computeForAreaInterior(geomA, tolerance); if (maxDistPtsInterior != null @@ -253,6 +350,7 @@ && distance(maxDistPtsInterior) > distance(maxDistPtsEdge)) { return maxDistPtsInterior; } } + //*/ return maxDistPtsEdge; } @@ -423,15 +521,23 @@ private Coordinate[] computeForAreaInterior(Geometry geomA, double tolerance) { /** * Optimization - skip if A interior cannot intersect B, - * and thus farther point must lie on A segment + * and thus farther point must lie on A boundary */ if (polygonalA.getEnvelopeInternal().disjoint(geomB.getEnvelopeInternal())) { return null; } + /** + * The LargestEmptyCircle computation is much slower than the boundary one, + * is quite unlikely to occur, + * and accuracy is probably less critical (or obvious). + * So improve performance by using a coarser distance tolerance. + */ + double lecTol = AREA_INTERIOR_PERFORMANCE_FACTOR * tolerance; + //TODO: add short-circuiting based on maxDistanceLimit? - Point centerPt = LargestEmptyCircle.getCenter(geomB, polygonalA, tolerance); + Point centerPt = LargestEmptyCircle.getCenter(geomB, polygonalA, lecTol); Coordinate ptA = centerPt.getCoordinate(); /** * If LEC centre is in B, the max distance is zero, so return null. diff --git a/modules/core/src/test/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistanceTest.java b/modules/core/src/test/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistanceTest.java index ab6c1880df..5971481558 100644 --- a/modules/core/src/test/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistanceTest.java +++ b/modules/core/src/test/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistanceTest.java @@ -114,7 +114,7 @@ public void testPolygonLineCrossingInteriorPoint() checkDistanceStartPtLen("POLYGON ((2 8, 8 2, 2 1, 2 8))", "LINESTRING (6 5, 4 7, 0 0, 9 1)", 0.001, - "LINESTRING (4.557 2.9937, 2.4 4.2)", 0.01); + "LINESTRING (4.555 2.989, 4.828 0.536)", 0.01); } public void testPolygonPolygon() From 4dd7927346f177f02e48cc38cad679fc219a0218 Mon Sep 17 00:00:00 2001 From: Martin Davis Date: Mon, 23 Feb 2026 13:14:14 -0800 Subject: [PATCH 31/43] Improve search termination criteria --- .../distance/DirectedHausdorffDistance.java | 27 ++++++++++++------- 1 file changed, 17 insertions(+), 10 deletions(-) diff --git a/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java b/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java index d8ea1344c4..ea55e5d78e 100644 --- a/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java +++ b/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java @@ -342,7 +342,6 @@ private Coordinate[] computeDistancePoints(Geometry geomA, double tolerance, dou /** * Polygonal query geometry may have an interior point as the farthest point. */ - //* if (geomA.getDimension() == Dimension.A) { Coordinate[] maxDistPtsInterior = computeForAreaInterior(geomA, tolerance); if (maxDistPtsInterior != null @@ -350,7 +349,6 @@ && distance(maxDistPtsInterior) > distance(maxDistPtsEdge)) { return maxDistPtsInterior; } } - //*/ return maxDistPtsEdge; } @@ -393,6 +391,9 @@ private Coordinate[] computeForEdges(Geometry geomA, double tolerance, double ma iter++; // get the segment with greatest distance bound DHDSegment segMaxBound = segQueue.remove(); + //System.out.println(segMaxBound); + //System.out.println(WKTWriter.toLineString(segMaxBound.getMaxDistPts())); + /* double maxDistBound = segMaxBound.maxDistanceBound(); double maxDist = segMaxBound.maxDistance(); @@ -403,14 +404,21 @@ private Coordinate[] computeForEdges(Geometry geomA, double tolerance, double ma } */ /** - * Save if segment is farther than most distant so far + * Save if segment point is farther than current farthest */ if (segMaxDist == null || segMaxBound.getMaxDistance() > segMaxDist.getMaxDistance()) { segMaxDist = segMaxBound; } /** - * If maxDistanceLimit is specified, short-circuit if: + * Stop searching if remaining items in queue must all be closer + * than the current maximum distance. + */ + if (segMaxBound.getMaxDistanceBound() <= segMaxDist.getMaxDistance()) { + break; + } + /** + * If maxDistanceLimit is specified, can stop searching if: * - if segment distance bound is less than distance limit, no other segment can be farther * - if a point of segment is farther than limit, isFulyWithin must be false */ @@ -419,7 +427,7 @@ private Coordinate[] computeForEdges(Geometry geomA, double tolerance, double ma ) { break; } - + /** * Check for equal or collinear segments. * If so, don't bisect the segment further. @@ -435,12 +443,11 @@ private Coordinate[] computeForEdges(Geometry geomA, double tolerance, double ma //System.out.println(segMaxDist); /** - * If segment is longer than tolerance - * and it might provide a better max distance point, - * bisect and keep searching + * If segment is longer than tolerance + * it might provide a better max distance point, + * so bisect and keep searching */ - if ((segMaxBound.getLength() > tolerance) - && segMaxBound.getMaxDistanceBound() > segMaxDist.getMaxDistance()) { + if ((segMaxBound.getLength() > tolerance)) { DHDSegment[] bisects = segMaxBound.bisect(distanceToB); addNonInterior(bisects[0], segQueue); addNonInterior(bisects[1], segQueue); From b8f8c4d0b70994658c1c6e67fdc1c6b1de6d1927 Mon Sep 17 00:00:00 2001 From: Martin Davis Date: Mon, 23 Feb 2026 17:51:56 -0800 Subject: [PATCH 32/43] Renames, remove hausdorff distance function with tolerance --- .../jtstest/function/DistanceFunctions.java | 6 +- .../distance/DirectedHausdorffDistance.java | 60 +++++++------------ .../DirectedHausdorffDistanceTest.java | 4 +- .../distance/HausdorffDistancePerfTest.java | 2 +- .../distance/HausdorffDistanceStressTest.java | 2 +- 5 files changed, 27 insertions(+), 47 deletions(-) diff --git a/modules/app/src/main/java/org/locationtech/jtstest/function/DistanceFunctions.java b/modules/app/src/main/java/org/locationtech/jtstest/function/DistanceFunctions.java index 71e6a177df..be0294d03e 100644 --- a/modules/app/src/main/java/org/locationtech/jtstest/function/DistanceFunctions.java +++ b/modules/app/src/main/java/org/locationtech/jtstest/function/DistanceFunctions.java @@ -97,11 +97,9 @@ public static Geometry directedHausdorffLine(Geometry a, Geometry b) } @Metadata(description="Hausdorff distance between A and B, up to tolerance") - public static Geometry hausdorffLine(Geometry a, Geometry b, - @Metadata(title="Distance tolerance") - double distTol) + public static Geometry hausdorffLine(Geometry a, Geometry b) { - return DirectedHausdorffDistance.hausdorffDistanceLine(a, b, distTol); + return DirectedHausdorffDistance.hausdorffDistanceLine(a, b); } //-------------------------------------------- diff --git a/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java b/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java index ea55e5d78e..bf54af88ca 100644 --- a/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java +++ b/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java @@ -69,7 +69,7 @@ *

* The class can be used in prepared mode. * Creating an instance on a target geometry caches indexes for that geometry. - * Then {@link #maximumDistancePoints(Geometry, double) + * Then {@link #farthestPoints(Geometry, double) * or {@link #isFullyWithinDistance(Geometry, double, double)} * can be called efficiently for multiple query geometries. *

@@ -98,7 +98,7 @@ public class DirectedHausdorffDistance { public static double distance(Geometry a, Geometry b, double tolerance) { DirectedHausdorffDistance hd = new DirectedHausdorffDistance(b); - return distance(hd.maximumDistancePoints(a, tolerance)); + return distance(hd.farthestPoints(a, tolerance)); } /** @@ -112,7 +112,7 @@ public static double distance(Geometry a, Geometry b, double tolerance) public static double distance(Geometry a, Geometry b) { DirectedHausdorffDistance hd = new DirectedHausdorffDistance(b); - return distance(hd.maximumDistancePoints(a)); + return distance(hd.farthestPoints(a)); } /** @@ -127,7 +127,7 @@ public static double distance(Geometry a, Geometry b) public static LineString distanceLine(Geometry a, Geometry b, double tolerance) { DirectedHausdorffDistance hd = new DirectedHausdorffDistance(b); - return a.getFactory().createLineString(hd.maximumDistancePoints(a, tolerance)); + return a.getFactory().createLineString(hd.farthestPoints(a, tolerance)); } /** @@ -142,33 +142,7 @@ public static LineString distanceLine(Geometry a, Geometry b, double tolerance) public static LineString distanceLine(Geometry a, Geometry b) { DirectedHausdorffDistance hd = new DirectedHausdorffDistance(b); - return a.getFactory().createLineString(hd.maximumDistancePoints(a)); - } - - /** - * Computes a pair of points which attain the symmetric Hausdorff distance - * between two geometries, up to a given distance accuracy. - * This the maximum of the two directed Hausdorff distances. - * - * @param a a geometry - * @param b a geometry - * @param tolerance the approximation distance tolerance - * @return a pair of points [ptA, ptB] demonstrating the Hausdorff distance - */ - public static LineString hausdorffDistanceLine(Geometry a, Geometry b, double tolerance) - { - DirectedHausdorffDistance hdAB = new DirectedHausdorffDistance(b); - Coordinate[] ptsAB = hdAB.maximumDistancePoints(a, tolerance); - DirectedHausdorffDistance hdBA = new DirectedHausdorffDistance(a); - Coordinate[] ptsBA = hdBA.maximumDistancePoints(b, tolerance); - - //-- return points in A-B order - Coordinate[] pts = ptsAB; - if (distance(ptsBA) > distance(ptsAB)) { - //-- reverse the BA points - pts = pair(ptsBA[1], ptsBA[0]); - } - return a.getFactory().createLineString(pts); + return a.getFactory().createLineString(hd.farthestPoints(a)); } /** @@ -183,9 +157,9 @@ public static LineString hausdorffDistanceLine(Geometry a, Geometry b, double to public static LineString hausdorffDistanceLine(Geometry a, Geometry b) { DirectedHausdorffDistance hdAB = new DirectedHausdorffDistance(b); - Coordinate[] ptsAB = hdAB.maximumDistancePoints(a); + Coordinate[] ptsAB = hdAB.farthestPoints(a); DirectedHausdorffDistance hdBA = new DirectedHausdorffDistance(a); - Coordinate[] ptsBA = hdBA.maximumDistancePoints(b); + Coordinate[] ptsBA = hdBA.farthestPoints(b); //-- return points in A-B order Coordinate[] pts = ptsAB; @@ -236,7 +210,7 @@ private static Coordinate[] pair(Coordinate p0, Coordinate p1) { /** * Heuristic automatic tolerance factor */ - private static final double AUTO_TOLERANCE_FACTOR = 1.0e3; + private static final double AUTO_TOLERANCE_FACTOR = 1.0e4; private static double computeTolerance(Geometry geom) { return geom.getEnvelopeInternal().getDiameter() / AUTO_TOLERANCE_FACTOR; @@ -311,20 +285,28 @@ private static boolean isBeyond(Envelope envA, Envelope envB, double maxDistance || envA.getMaxY() > envB.getMaxY() + maxDistance; } - public Coordinate[] maximumDistancePoints(Geometry geomA) { + /** + * Computes a pair of points which attain the directed Hausdorff distance + * of a query geometry A from the target B. + * + * @param geomA the query geometry + * @return a pair of points [ptA, ptB] attaining the distance + */ + public Coordinate[] farthestPoints(Geometry geomA) { double tolerance = computeTolerance(geomA); - return computeDistancePoints(geomA, tolerance, -1.0); + return farthestPoints(geomA, tolerance); } /** * Computes a pair of points which attain the directed Hausdorff distance - * of a query geometry A from the target B. + * of a query geometry A from the target B, + * up to a given distance accuracy. * * @param geomA the query geometry * @param tolerance the approximation distance tolerance - * @return a pair of points [ptA, ptB] demonstrating the distance + * @return a pair of points [ptA, ptB] attaining the distance */ - public Coordinate[] maximumDistancePoints(Geometry geomA, double tolerance) { + public Coordinate[] farthestPoints(Geometry geomA, double tolerance) { return computeDistancePoints(geomA, tolerance, -1.0); } diff --git a/modules/core/src/test/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistanceTest.java b/modules/core/src/test/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistanceTest.java index 5971481558..d8ff2556eb 100644 --- a/modules/core/src/test/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistanceTest.java +++ b/modules/core/src/test/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistanceTest.java @@ -158,7 +158,7 @@ public void testNonVertexResult() String wkt1 = "LINESTRING (1 1, 5 10, 9 1)"; String wkt2 = "LINESTRING (0 10, 0 0, 10 0)"; - checkHD(wkt1, wkt2, 0.01, "LINESTRING (6.5390625 6.537109375, 6.5390625 0)"); + checkHD(wkt1, wkt2, 0.01, "LINESTRING (6.53857421875 6.5382080078125, 6.53857421875 0)"); checkDistance(wkt1, wkt2, 0.01, "LINESTRING (6.5390625 6.537109375, 6.5390625 0)"); } @@ -266,7 +266,7 @@ private void checkHD(String wkt1, String wkt2, double tolerance, String wktExpec Geometry g1 = read(wkt1); Geometry g2 = read(wkt2); - Geometry result = DirectedHausdorffDistance.hausdorffDistanceLine(g1, g2, tolerance); + Geometry result = DirectedHausdorffDistance.hausdorffDistanceLine(g1, g2); Geometry expected = read(wktExpected); checkEqualExact(expected, result, TOLERANCE); diff --git a/modules/core/src/test/java/test/jts/perf/operation/distance/HausdorffDistancePerfTest.java b/modules/core/src/test/java/test/jts/perf/operation/distance/HausdorffDistancePerfTest.java index c6582afb13..9b4c113d5f 100644 --- a/modules/core/src/test/java/test/jts/perf/operation/distance/HausdorffDistancePerfTest.java +++ b/modules/core/src/test/java/test/jts/perf/operation/distance/HausdorffDistancePerfTest.java @@ -82,7 +82,7 @@ public void runFullyWithin() { private void checkHausdorff(Geometry g, Geometry b) { double distDHD = DiscreteHausdorffDistance.distance(g, b, 0.01); - double distHD = DirectedHausdorffDistance.hausdorffDistanceLine(g, b, 0.01).getLength(); + double distHD = DirectedHausdorffDistance.hausdorffDistanceLine(g, b).getLength(); //-- performance testing only //double distDHD = distHD; diff --git a/modules/core/src/test/java/test/jts/perf/operation/distance/HausdorffDistanceStressTest.java b/modules/core/src/test/java/test/jts/perf/operation/distance/HausdorffDistanceStressTest.java index a1513ffe9a..9fa4d1b98e 100644 --- a/modules/core/src/test/java/test/jts/perf/operation/distance/HausdorffDistanceStressTest.java +++ b/modules/core/src/test/java/test/jts/perf/operation/distance/HausdorffDistanceStressTest.java @@ -56,7 +56,7 @@ private void run() { private void checkHausdorff(Geometry g, Geometry b) { double distDHD = DiscreteHausdorffDistance.distance(g, b, 0.01); - double distHD = DirectedHausdorffDistance.hausdorffDistanceLine(g, b, 0.01).getLength(); + double distHD = DirectedHausdorffDistance.hausdorffDistanceLine(g, b).getLength(); //-- performance testing only //double distDHD = distHD; From e6d121000376d0694437b0fac46bf417edcfa252 Mon Sep 17 00:00:00 2001 From: Martin Davis Date: Mon, 23 Feb 2026 20:33:36 -0800 Subject: [PATCH 33/43] Remove tolerance from Hausdorff test --- .../distance/DirectedHausdorffDistanceTest.java | 15 +++++---------- 1 file changed, 5 insertions(+), 10 deletions(-) diff --git a/modules/core/src/test/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistanceTest.java b/modules/core/src/test/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistanceTest.java index d8ff2556eb..86a08d3e5f 100644 --- a/modules/core/src/test/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistanceTest.java +++ b/modules/core/src/test/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistanceTest.java @@ -29,7 +29,6 @@ public static void main(String args[]) { public void testPointPoint() { checkHD("POINT (0 0)", "POINT (1 1)", - 0.01, "LINESTRING (0 0, 1 1)"); } @@ -39,7 +38,7 @@ public void testPointsPoints() String b = "MULTIPOINT ((0.1 0), (1 0), (2 0), (3 0), (4 0), (5 0))"; checkDistance(a, b, 0.01, "LINESTRING (6 6, 5 0)"); checkDistance(b, a, 0.01, "LINESTRING (5 0, 2 3)"); - checkHD(a, b, 0.01, "LINESTRING (6 6, 5 0)"); + checkHD(a, b, "LINESTRING (6 6, 5 0)"); } public void testPointPolygonInterior() @@ -59,21 +58,18 @@ public void testPointsPolygon() public void testLineSegments() { checkHD("LINESTRING (0 0, 2 0)", "LINESTRING (0 0, 2 1)", - 0.01, "LINESTRING (2 0, 2 1)"); } public void testLineSegments2() { checkHD("LINESTRING (0 0, 2 0)", "LINESTRING (0 1, 1 2, 2 1)", - 0.01, "LINESTRING (1 0, 1 2)"); } public void testLinePoints() { checkHD("LINESTRING (0 0, 2 0)", "MULTIPOINT (0 2, 1 0, 2 1)", - 0.01, "LINESTRING (0 0, 0 2)"); } @@ -89,7 +85,6 @@ public void testLinesPolygon() { checkHD("MULTILINESTRING ((1 1, 2 7), (7 1, 9 9))", "POLYGON ((3 7, 6 7, 6 4, 3 4, 3 7))", - 0.01, "LINESTRING (9 9, 6 7)"); } @@ -98,7 +93,7 @@ public void testLinesPolygon2() String a = "MULTILINESTRING ((2 3, 2 7), (9 1, 9 8, 4 9))"; String b = "POLYGON ((3 7, 6 8, 8 2, 3 4, 3 7))"; checkDistance(a, b, 0.01, "LINESTRING (9 8, 6.3 7.1)"); - checkHD(a, b, 0.01, "LINESTRING (2 3, 5.5 3)"); + checkHD(a, b, "LINESTRING (2 3, 5.5 3)"); } public void testPolygonLineCrossingBoundaryResult() @@ -123,7 +118,7 @@ public void testPolygonPolygon() String b = "POLYGON ((1 19, 5 12, 5 3, 14 10, 11 19, 19 19, 20 0, 1 1, 1 19))"; checkDistance(b, a, 0.01, "LINESTRING (20 0, 17 3)"); checkDistance(a, b, 0.01, "LINESTRING (6.6796875 18, 11 19)"); - checkHD(a, b, 0.01, "LINESTRING (6.6796875 18, 11 19)"); + checkHD(a, b, "LINESTRING (6.6796875 18, 11 19)"); } public void testPolygonPolygonHolesNested() @@ -158,7 +153,7 @@ public void testNonVertexResult() String wkt1 = "LINESTRING (1 1, 5 10, 9 1)"; String wkt2 = "LINESTRING (0 10, 0 0, 10 0)"; - checkHD(wkt1, wkt2, 0.01, "LINESTRING (6.53857421875 6.5382080078125, 6.53857421875 0)"); + checkHD(wkt1, wkt2, "LINESTRING (6.53857421875 6.5382080078125, 6.53857421875 0)"); checkDistance(wkt1, wkt2, 0.01, "LINESTRING (6.5390625 6.537109375, 6.5390625 0)"); } @@ -261,7 +256,7 @@ public void testFullyWithinDistancePolygonsNestedWithHole() private static final double TOLERANCE = 0.001; - private void checkHD(String wkt1, String wkt2, double tolerance, String wktExpected) + private void checkHD(String wkt1, String wkt2, String wktExpected) { Geometry g1 = read(wkt1); Geometry g2 = read(wkt2); From 8480ff4165a23e8ad528e000148ead7a64135eb9 Mon Sep 17 00:00:00 2001 From: Martin Davis Date: Mon, 23 Feb 2026 20:36:47 -0800 Subject: [PATCH 34/43] Remove tolerance from fullyWithin test --- .../DirectedHausdorffDistanceTest.java | 68 +++++++++---------- 1 file changed, 31 insertions(+), 37 deletions(-) diff --git a/modules/core/src/test/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistanceTest.java b/modules/core/src/test/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistanceTest.java index 86a08d3e5f..f905fac90f 100644 --- a/modules/core/src/test/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistanceTest.java +++ b/modules/core/src/test/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistanceTest.java @@ -28,7 +28,7 @@ public static void main(String args[]) { public void testPointPoint() { - checkHD("POINT (0 0)", "POINT (1 1)", + checkHausdorff("POINT (0 0)", "POINT (1 1)", "LINESTRING (0 0, 1 1)"); } @@ -38,7 +38,7 @@ public void testPointsPoints() String b = "MULTIPOINT ((0.1 0), (1 0), (2 0), (3 0), (4 0), (5 0))"; checkDistance(a, b, 0.01, "LINESTRING (6 6, 5 0)"); checkDistance(b, a, 0.01, "LINESTRING (5 0, 2 3)"); - checkHD(a, b, "LINESTRING (6 6, 5 0)"); + checkHausdorff(a, b, "LINESTRING (6 6, 5 0)"); } public void testPointPolygonInterior() @@ -57,19 +57,19 @@ public void testPointsPolygon() public void testLineSegments() { - checkHD("LINESTRING (0 0, 2 0)", "LINESTRING (0 0, 2 1)", + checkHausdorff("LINESTRING (0 0, 2 0)", "LINESTRING (0 0, 2 1)", "LINESTRING (2 0, 2 1)"); } public void testLineSegments2() { - checkHD("LINESTRING (0 0, 2 0)", "LINESTRING (0 1, 1 2, 2 1)", + checkHausdorff("LINESTRING (0 0, 2 0)", "LINESTRING (0 1, 1 2, 2 1)", "LINESTRING (1 0, 1 2)"); } public void testLinePoints() { - checkHD("LINESTRING (0 0, 2 0)", "MULTIPOINT (0 2, 1 0, 2 1)", + checkHausdorff("LINESTRING (0 0, 2 0)", "MULTIPOINT (0 2, 1 0, 2 1)", "LINESTRING (0 0, 0 2)"); } @@ -83,7 +83,7 @@ public void testLinesTopoEqual() public void testLinesPolygon() { - checkHD("MULTILINESTRING ((1 1, 2 7), (7 1, 9 9))", + checkHausdorff("MULTILINESTRING ((1 1, 2 7), (7 1, 9 9))", "POLYGON ((3 7, 6 7, 6 4, 3 4, 3 7))", "LINESTRING (9 9, 6 7)"); } @@ -93,7 +93,7 @@ public void testLinesPolygon2() String a = "MULTILINESTRING ((2 3, 2 7), (9 1, 9 8, 4 9))"; String b = "POLYGON ((3 7, 6 8, 8 2, 3 4, 3 7))"; checkDistance(a, b, 0.01, "LINESTRING (9 8, 6.3 7.1)"); - checkHD(a, b, "LINESTRING (2 3, 5.5 3)"); + checkHausdorff(a, b, "LINESTRING (2 3, 5.5 3)"); } public void testPolygonLineCrossingBoundaryResult() @@ -118,7 +118,7 @@ public void testPolygonPolygon() String b = "POLYGON ((1 19, 5 12, 5 3, 14 10, 11 19, 19 19, 20 0, 1 1, 1 19))"; checkDistance(b, a, 0.01, "LINESTRING (20 0, 17 3)"); checkDistance(a, b, 0.01, "LINESTRING (6.6796875 18, 11 19)"); - checkHD(a, b, "LINESTRING (6.6796875 18, 11 19)"); + checkHausdorff(a, b, "LINESTRING (6.6796875 18, 11 19)"); } public void testPolygonPolygonHolesNested() @@ -153,7 +153,7 @@ public void testNonVertexResult() String wkt1 = "LINESTRING (1 1, 5 10, 9 1)"; String wkt2 = "LINESTRING (0 10, 0 0, 10 0)"; - checkHD(wkt1, wkt2, "LINESTRING (6.53857421875 6.5382080078125, 6.53857421875 0)"); + checkHausdorff(wkt1, wkt2, "LINESTRING (6.53857421875 6.5382080078125, 6.53857421875 0)"); checkDistance(wkt1, wkt2, 0.01, "LINESTRING (6.5390625 6.537109375, 6.5390625 0)"); } @@ -201,62 +201,62 @@ public void testFullyWithinDistancePoints() { String a = "MULTIPOINT ((1 9), (9 1))"; String b = "MULTIPOINT ((1 1), (9 9))"; - checkFullyWithinDistance(a, b, 1, 0.01, false); - checkFullyWithinDistance(a, b, 8.1, 0.01, true); + checkFullyWithinDistance(a, b, 1, false); + checkFullyWithinDistance(a, b, 8.1, true); } public void testFullyWithinDistanceDisconnectedLines() { String a = "MULTILINESTRING ((1 9, 2 9), (8 1, 9 1))"; String b = "LINESTRING (9 9, 1 1)"; - checkFullyWithinDistance(a, b, 1, 0.01, false); - checkFullyWithinDistance(a, b, 6, 0.01, true); - checkFullyWithinDistance(b, a, 1, 0.01, false); - checkFullyWithinDistance(b, a, 7.1, 0.01, true); + checkFullyWithinDistance(a, b, 1, false); + checkFullyWithinDistance(a, b, 6, true); + checkFullyWithinDistance(b, a, 1, false); + checkFullyWithinDistance(b, a, 7.1, true); } public void testFullyWithinDistanceDisconnectedPolygons() { String a = "MULTIPOLYGON (((1 9, 2 9, 2 8, 1 8, 1 9)), ((8 2, 9 2, 9 1, 8 1, 8 2)))"; String b = "POLYGON ((1 2, 9 9, 2 1, 1 2))"; - checkFullyWithinDistance(a, b, 1, 0.01, false); - checkFullyWithinDistance(a, b, 5.3, 0.01, true); - checkFullyWithinDistance(b, a, 1, 0.01, false); - checkFullyWithinDistance(b, a, 7.1, 0.01, true); + checkFullyWithinDistance(a, b, 1, false); + checkFullyWithinDistance(a, b, 5.3, true); + checkFullyWithinDistance(b, a, 1, false); + checkFullyWithinDistance(b, a, 7.1, true); } public void testFullyWithinDistanceLines() { String a = "MULTILINESTRING ((1 1, 3 3), (7 7, 9 9))"; String b = "MULTILINESTRING ((1 9, 1 5), (6 4, 8 2))"; - checkFullyWithinDistance(a, b, 1, 0.01, false); - checkFullyWithinDistance(a, b, 4, 0.01, false); - checkFullyWithinDistance(a, b, 6, 0.01, true); + checkFullyWithinDistance(a, b, 1, false); + checkFullyWithinDistance(a, b, 4, false); + checkFullyWithinDistance(a, b, 6, true); } public void testFullyWithinDistancePolygons() { String a = "POLYGON ((1 4, 4 4, 4 1, 1 1, 1 4))"; String b = "POLYGON ((10 10, 10 15, 15 15, 15 10, 10 10))"; - checkFullyWithinDistance(a, b, 5, 0.01, false); - checkFullyWithinDistance(a, b, 10, 0.01, false); - checkFullyWithinDistance(a, b, 20, 0.01, true); + checkFullyWithinDistance(a, b, 5, false); + checkFullyWithinDistance(a, b, 10, false); + checkFullyWithinDistance(a, b, 20, true); } public void testFullyWithinDistancePolygonsNestedWithHole() { String a = "POLYGON ((2 8, 8 8, 8 2, 2 2, 2 8))"; String b = "POLYGON ((1 9, 9 9, 9 1, 1 1, 1 9), (3 7, 7 7, 7 3, 3 3, 3 7))"; - checkFullyWithinDistance(a, b, 1, 0.01, false); - checkFullyWithinDistance(a, b, 2, 0.01, true); - checkFullyWithinDistance(a, b, 3, 0.01, true); + checkFullyWithinDistance(a, b, 1, false); + checkFullyWithinDistance(a, b, 2, true); + checkFullyWithinDistance(a, b, 3, true); } //====================================================================== private static final double TOLERANCE = 0.001; - private void checkHD(String wkt1, String wkt2, String wktExpected) + private void checkHausdorff(String wkt1, String wkt2, String wktExpected) { Geometry g1 = read(wkt1); Geometry g2 = read(wkt2); @@ -264,12 +264,6 @@ private void checkHD(String wkt1, String wkt2, String wktExpected) Geometry result = DirectedHausdorffDistance.hausdorffDistanceLine(g1, g2); Geometry expected = read(wktExpected); checkEqualExact(expected, result, TOLERANCE); - - /* - double resultDistance = DiscreteHausdorffDistance.distance(g1, g2); - double expectedDistance = expected.getLength(); - assertEquals(expectedDistance, resultDistance, TOLERANCE); - */ } private void checkDistance(String wkt1, String wkt2, double tolerance, String wktExpected) { @@ -310,11 +304,11 @@ private void checkDistance(String wkt1, String wkt2, double tolerance, } - private void checkFullyWithinDistance(String a, String b, double distance, double tolerance, boolean expected) { + private void checkFullyWithinDistance(String a, String b, double distance, boolean expected) { Geometry g1 = read(a); Geometry g2 = read(b); - boolean result = DirectedHausdorffDistance.isFullyWithinDistance(g1, g2, distance, tolerance); + boolean result = DirectedHausdorffDistance.isFullyWithinDistance(g1, g2, distance); assertEquals(expected, result); } From 811e3205d3a242e90aa669a822f7c2b72a413985 Mon Sep 17 00:00:00 2001 From: Martin Davis Date: Mon, 23 Feb 2026 21:57:53 -0800 Subject: [PATCH 35/43] Remove tolerance from tests --- .../DirectedHausdorffDistanceTest.java | 62 ++++++++++--------- 1 file changed, 33 insertions(+), 29 deletions(-) diff --git a/modules/core/src/test/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistanceTest.java b/modules/core/src/test/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistanceTest.java index f905fac90f..9b94c525a5 100644 --- a/modules/core/src/test/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistanceTest.java +++ b/modules/core/src/test/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistanceTest.java @@ -36,22 +36,20 @@ public void testPointsPoints() { String a = "MULTIPOINT ((0 1), (2 3), (4 5), (6 6))"; String b = "MULTIPOINT ((0.1 0), (1 0), (2 0), (3 0), (4 0), (5 0))"; - checkDistance(a, b, 0.01, "LINESTRING (6 6, 5 0)"); - checkDistance(b, a, 0.01, "LINESTRING (5 0, 2 3)"); + checkDistance(a, b, "LINESTRING (6 6, 5 0)"); + checkDistance(b, a, "LINESTRING (5 0, 2 3)"); checkHausdorff(a, b, "LINESTRING (6 6, 5 0)"); } public void testPointPolygonInterior() { checkDistance("POINT (3 4)", "POLYGON ((1 9, 9 9, 9 1, 1 1, 1 9))", - 0.01, 0); } public void testPointsPolygon() { checkDistance("MULTIPOINT ((4 3), (2 8), (8 5))", "POLYGON ((6 9, 6 4, 9 1, 1 1, 6 9))", - 0.01, "LINESTRING (2 8, 4.426966292134832 6.48314606741573)"); } @@ -77,7 +75,6 @@ public void testLinesTopoEqual() { checkDistance("MULTILINESTRING ((10 10, 10 90, 40 30), (40 30, 60 80, 90 30, 40 10))", "LINESTRING (10 10, 10 90, 40 30, 60 80, 90 30, 40 10)", - 0.01, 0.0); } @@ -92,7 +89,7 @@ public void testLinesPolygon2() { String a = "MULTILINESTRING ((2 3, 2 7), (9 1, 9 8, 4 9))"; String b = "POLYGON ((3 7, 6 8, 8 2, 3 4, 3 7))"; - checkDistance(a, b, 0.01, "LINESTRING (9 8, 6.3 7.1)"); + checkDistance(a, b, "LINESTRING (9 8, 6.3 7.1)"); checkHausdorff(a, b, "LINESTRING (2 3, 5.5 3)"); } @@ -100,7 +97,6 @@ public void testPolygonLineCrossingBoundaryResult() { checkDistance("POLYGON ((2 8, 8 2, 2 1, 2 8))", "LINESTRING (6 5, 4 7, 0 0, 8 4)", - 0.001, "LINESTRING (2 8, 3.9384615384615387 6.892307692307693)"); } @@ -108,7 +104,6 @@ public void testPolygonLineCrossingInteriorPoint() { checkDistanceStartPtLen("POLYGON ((2 8, 8 2, 2 1, 2 8))", "LINESTRING (6 5, 4 7, 0 0, 9 1)", - 0.001, "LINESTRING (4.555 2.989, 4.828 0.536)", 0.01); } @@ -116,8 +111,8 @@ public void testPolygonPolygon() { String a = "POLYGON ((2 18, 18 18, 17 3, 2 2, 2 18))"; String b = "POLYGON ((1 19, 5 12, 5 3, 14 10, 11 19, 19 19, 20 0, 1 1, 1 19))"; - checkDistance(b, a, 0.01, "LINESTRING (20 0, 17 3)"); - checkDistance(a, b, 0.01, "LINESTRING (6.6796875 18, 11 19)"); + checkDistance(b, a, "LINESTRING (20 0, 17 3)"); + checkDistance(a, b, "LINESTRING (6.6796875 18, 11 19)"); checkHausdorff(a, b, "LINESTRING (6.6796875 18, 11 19)"); } @@ -126,16 +121,16 @@ public void testPolygonPolygonHolesNested() // B is contained in A String a = "POLYGON ((1 19, 19 19, 19 1, 1 1, 1 19), (6 8, 11 14, 15 7, 6 8))"; String b = "POLYGON ((2 18, 18 18, 18 2, 2 2, 2 18), (10 17, 3 7, 17 5, 10 17))"; - checkDistance(a, b, 0.01, "LINESTRING (9.8134765625 12.576171875, 7.8603515625 13.943359375)"); - checkDistance(b, a, 0.01, 0.0); + checkDistance(a, b, "LINESTRING (9.817138671875 12.58056640625, 7.863620425230705 13.948029178901006)"); + checkDistance(b, a, 0.0); } public void testMultiPolygons() { String a = "MULTIPOLYGON (((1 1, 1 10, 5 1, 1 1)), ((4 17, 9 15, 9 6, 4 17)))"; String b = "MULTIPOLYGON (((1 12, 4 13, 8 10, 1 12)), ((3 8, 7 7, 6 2, 3 8)))"; - checkDistance(a, b, 0.01, "LINESTRING (1 1, 5.4 3.2)"); - checkDistanceStartPtLen(b, a, 0.001, + checkDistance(a, b, "LINESTRING (1 1, 5.4 3.2)"); + checkDistanceStartPtLen(b, a, "LINESTRING (2.669921875 12.556640625, 5.446115154109589 13.818546660958905)", 0.01); } @@ -145,7 +140,7 @@ public void testLinePolygonCrossing() throws Exception { String wkt1 = "LINESTRING (2 5, 5 10, 6 4)"; String wkt2 = "POLYGON ((1 9, 9 9, 9 1, 1 1, 1 9))"; - checkDistance(wkt1, wkt2, 0.01, "LINESTRING (5 10, 5 9)"); + checkDistance(wkt1, wkt2, "LINESTRING (5 10, 5 9)"); } public void testNonVertexResult() @@ -154,23 +149,23 @@ public void testNonVertexResult() String wkt2 = "LINESTRING (0 10, 0 0, 10 0)"; checkHausdorff(wkt1, wkt2, "LINESTRING (6.53857421875 6.5382080078125, 6.53857421875 0)"); - checkDistance(wkt1, wkt2, 0.01, "LINESTRING (6.5390625 6.537109375, 6.5390625 0)"); + checkDistance(wkt1, wkt2, "LINESTRING (6.53857421875 6.5382080078125, 6.53857421875 0)"); } public void testDirectedLines() { String wkt1 = "LINESTRING (1 6, 3 5, 1 4)"; String wkt2 = "LINESTRING (1 10, 9 5, 1 2)"; - checkDistance(wkt1, wkt2, 0.01, "LINESTRING (1 6, 2.797752808988764 8.876404494382022)"); - checkDistance(wkt2, wkt1, 0.01, "LINESTRING (9 5, 3 5)"); + checkDistance(wkt1, wkt2, "LINESTRING (1 6, 2.797752808988764 8.876404494382022)"); + checkDistance(wkt2, wkt1, "LINESTRING (9 5, 3 5)"); } public void testDirectedLines2() { String wkt1 = "LINESTRING (1 6, 3 5, 1 4)"; String wkt2 = "LINESTRING (1 3, 1 9, 9 5, 1 1)"; - checkDistance(wkt1, wkt2, 0.01, "LINESTRING (3 5, 1 5)"); - checkDistance(wkt2, wkt1, 0.01, "LINESTRING (9 5, 3 5)"); + checkDistance(wkt1, wkt2, "LINESTRING (3 5, 1 5)"); + checkDistance(wkt2, wkt1, "LINESTRING (9 5, 3 5)"); } /** @@ -191,7 +186,7 @@ public void testInteriorSegmentsSameExterior() { String a = "POLYGON ((1 9, 3 9, 4 5, 5.05 9, 9 9, 9 1, 1 1, 1 9))"; String b = "POLYGON ((1 9, 9 9, 9 1, 1 1, 1 9))"; - checkDistance(a, b, 1.0, 0.0); + checkDistance(a, b, 0.0); } //----------------------------------------------------- @@ -266,21 +261,21 @@ private void checkHausdorff(String wkt1, String wkt2, String wktExpected) checkEqualExact(expected, result, TOLERANCE); } - private void checkDistance(String wkt1, String wkt2, double tolerance, String wktExpected) { + private void checkDistance(String wkt1, String wkt2, String wktExpected) { Geometry g1 = read(wkt1); Geometry g2 = read(wkt2); - Geometry result = DirectedHausdorffDistance.distanceLine(g1, g2, tolerance); + Geometry result = DirectedHausdorffDistance.distanceLine(g1, g2); Geometry expected = read(wktExpected); - checkEqualExact(expected, result, 100 * tolerance); + checkEqualExact(expected, result, TOLERANCE); } - - private void checkDistanceStartPtLen(String wkt1, String wkt2, double tolerance, + + private void checkDistanceStartPtLen(String wkt1, String wkt2, String wktExpected, double resultTolerance) { Geometry g1 = read(wkt1); Geometry g2 = read(wkt2); - Geometry result = DirectedHausdorffDistance.distanceLine(g1, g2, tolerance); + Geometry result = DirectedHausdorffDistance.distanceLine(g1, g2); Geometry expected = read(wktExpected); Coordinate resultPt = result.getCoordinates()[0]; @@ -303,7 +298,17 @@ private void checkDistance(String wkt1, String wkt2, double tolerance, assertEquals(expectedDistance, distResult, TOLERANCE); } - + private void checkDistance(String wkt1, String wkt2, + double expectedDistance) { + Geometry g1 = read(wkt1); + Geometry g2 = read(wkt2); + + Geometry result = DirectedHausdorffDistance.distanceLine(g1, g2); + + double distResult = result.getLength(); + assertEquals(expectedDistance, distResult, TOLERANCE); + } + private void checkFullyWithinDistance(String a, String b, double distance, boolean expected) { Geometry g1 = read(a); Geometry g2 = read(b); @@ -312,6 +317,5 @@ private void checkFullyWithinDistance(String a, String b, double distance, boole assertEquals(expected, result); } - } From 5206b498ccb610ebfc29078061a395ad070b846a Mon Sep 17 00:00:00 2001 From: Martin Davis Date: Mon, 23 Feb 2026 22:19:01 -0800 Subject: [PATCH 36/43] return point pairs instead of lines --- .../jtstest/function/DistanceFunctions.java | 9 +- .../distance/DirectedHausdorffDistance.java | 82 +++++++++++++------ .../DirectedHausdorffDistanceTest.java | 18 ++-- .../distance/HausdorffDistancePerfTest.java | 11 ++- .../distance/HausdorffDistanceStressTest.java | 16 ++-- 5 files changed, 82 insertions(+), 54 deletions(-) diff --git a/modules/app/src/main/java/org/locationtech/jtstest/function/DistanceFunctions.java b/modules/app/src/main/java/org/locationtech/jtstest/function/DistanceFunctions.java index be0294d03e..10d7aa226b 100644 --- a/modules/app/src/main/java/org/locationtech/jtstest/function/DistanceFunctions.java +++ b/modules/app/src/main/java/org/locationtech/jtstest/function/DistanceFunctions.java @@ -87,19 +87,22 @@ public static Geometry directedHausdorffLineTol(Geometry a, Geometry b, @Metadata(title="Distance tolerance") double distTol) { - return DirectedHausdorffDistance.distanceLine(a, b, distTol); + Coordinate[] pts = DirectedHausdorffDistance.distancePoints(a, b, distTol); + return a.getFactory().createLineString(pts); } @Metadata(description="Directed Hausdorff distance line from A to B") public static Geometry directedHausdorffLine(Geometry a, Geometry b) { - return DirectedHausdorffDistance.distanceLine(a, b); + Coordinate[] pts = DirectedHausdorffDistance.distancePoints(a, b); + return a.getFactory().createLineString(pts); } @Metadata(description="Hausdorff distance between A and B, up to tolerance") public static Geometry hausdorffLine(Geometry a, Geometry b) { - return DirectedHausdorffDistance.hausdorffDistanceLine(a, b); + Coordinate[] pts = DirectedHausdorffDistance.hausdorffDistancePoints(a, b); + return a.getFactory().createLineString(pts); } //-------------------------------------------- diff --git a/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java b/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java index bf54af88ca..cca9f4a5a4 100644 --- a/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java +++ b/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java @@ -84,65 +84,65 @@ * */ public class DirectedHausdorffDistance { - + /** * Computes the directed Hausdorff distance - * of a query geometry A from a target one B, - * up to a given distance accuracy. + * of a query geometry A from a target one B. * * @param a the query geometry * @param b the target geometry - * @param tolerance the accuracy distance tolerance * @return the directed Hausdorff distance */ - public static double distance(Geometry a, Geometry b, double tolerance) + public static double distance(Geometry a, Geometry b) { DirectedHausdorffDistance hd = new DirectedHausdorffDistance(b); - return distance(hd.farthestPoints(a, tolerance)); + return distance(hd.farthestPoints(a)); } /** * Computes the directed Hausdorff distance - * of a query geometry A from a target one B. + * of a query geometry A from a target one B, + * up to a given distance accuracy. * * @param a the query geometry * @param b the target geometry + * @param tolerance the accuracy distance tolerance * @return the directed Hausdorff distance */ - public static double distance(Geometry a, Geometry b) + public static double distance(Geometry a, Geometry b, double tolerance) { DirectedHausdorffDistance hd = new DirectedHausdorffDistance(b); - return distance(hd.farthestPoints(a)); + return distance(hd.farthestPoints(a, tolerance)); } - + /** * Computes a line containing a pair of points which attain the directed Hausdorff distance - * of a query geometry A from a target one B, up to a given distance accuracy. + * of a query geometry A from a target one B. * * @param a the query geometry * @param b the target geometry * @param tolerance the accuracy distance tolerance * @return a pair of points [ptA, ptB] demonstrating the distance */ - public static LineString distanceLine(Geometry a, Geometry b, double tolerance) + public static Coordinate[] distancePoints(Geometry a, Geometry b) { - DirectedHausdorffDistance hd = new DirectedHausdorffDistance(b); - return a.getFactory().createLineString(hd.farthestPoints(a, tolerance)); + DirectedHausdorffDistance dhd = new DirectedHausdorffDistance(b); + return dhd.farthestPoints(a); } - + /** * Computes a line containing a pair of points which attain the directed Hausdorff distance - * of a query geometry A from a target one B. + * of a query geometry A from a target one B, up to a given distance accuracy. * * @param a the query geometry * @param b the target geometry * @param tolerance the accuracy distance tolerance * @return a pair of points [ptA, ptB] demonstrating the distance */ - public static LineString distanceLine(Geometry a, Geometry b) + public static Coordinate[] distancePoints(Geometry a, Geometry b, double tolerance) { - DirectedHausdorffDistance hd = new DirectedHausdorffDistance(b); - return a.getFactory().createLineString(hd.farthestPoints(a)); + DirectedHausdorffDistance dhd = new DirectedHausdorffDistance(b); + return dhd.farthestPoints(a, tolerance); } /** @@ -154,7 +154,7 @@ public static LineString distanceLine(Geometry a, Geometry b) * @param b a geometry * @return a pair of points [ptA, ptB] demonstrating the Hausdorff distance */ - public static LineString hausdorffDistanceLine(Geometry a, Geometry b) + public static Coordinate[] hausdorffDistancePoints(Geometry a, Geometry b) { DirectedHausdorffDistance hdAB = new DirectedHausdorffDistance(b); Coordinate[] ptsAB = hdAB.farthestPoints(a); @@ -167,7 +167,20 @@ public static LineString hausdorffDistanceLine(Geometry a, Geometry b) //-- reverse the BA points pts = pair(ptsBA[1], ptsBA[0]); } - return a.getFactory().createLineString(pts); + return pts; + } + + /** + * Computes the symmetric Hausdorff distance between two geometries. + * This the maximum of the two directed Hausdorff distances. + * + * @param a a geometry + * @param b a geometry + * @return the Hausdorff distance + */ + public static double hausdorffDistance(Geometry a, Geometry b) + { + return distance(hausdorffDistancePoints(a, b)); } /** @@ -179,19 +192,31 @@ public static LineString hausdorffDistanceLine(Geometry a, Geometry b) * @param a the query geometry * @param b the target geometry * @param maxDistance the distance limit - * @param tolerance the approximation distance tolerance * @return true if the query geometry lies fully within the distance of the target */ - public static boolean isFullyWithinDistance(Geometry a, Geometry b, double maxDistance, double tolerance) + public static boolean isFullyWithinDistance(Geometry a, Geometry b, double maxDistance) { DirectedHausdorffDistance hd = new DirectedHausdorffDistance(b); - return hd.isFullyWithinDistance(a, maxDistance, tolerance); + return hd.isFullyWithinDistance(a, maxDistance); } - public static boolean isFullyWithinDistance(Geometry a, Geometry b, double maxDistance) + /** + * Computes whether a query geometry lies fully within a give distance of a target geometry, + * up to a given distance accuracy. + * Equivalently, detects whether any point of the query geometry is farther + * from the target than the specified distance. + * This is the case if DHD(A, B) > maxDistance. + * + * @param a the query geometry + * @param b the target geometry + * @param maxDistance the distance limit + * @param tolerance the accuracy distance tolerance + * @return true if the query geometry lies fully within the distance of the target + */ + public static boolean isFullyWithinDistance(Geometry a, Geometry b, double maxDistance, double tolerance) { DirectedHausdorffDistance hd = new DirectedHausdorffDistance(b); - return hd.isFullyWithinDistance(a, maxDistance); + return hd.isFullyWithinDistance(a, maxDistance, tolerance); } private static double distance(Coordinate[] pts) { @@ -240,7 +265,10 @@ public DirectedHausdorffDistance(Geometry b) { * @return true if the query geometry lies fully within the distance of the target */ public boolean isFullyWithinDistance(Geometry a, double maxDistance) { - return isFullyWithinDistance(a, maxDistance, computeTolerance(a)); + //TODO: should the tolerance be computed as a fraction of the maxDistance? + //return isFullyWithinDistance(a, maxDistance, computeTolerance(a)); + double tolerance = maxDistance / AUTO_TOLERANCE_FACTOR; + return isFullyWithinDistance(a, maxDistance, tolerance); } /** diff --git a/modules/core/src/test/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistanceTest.java b/modules/core/src/test/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistanceTest.java index 9b94c525a5..fa77d17670 100644 --- a/modules/core/src/test/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistanceTest.java +++ b/modules/core/src/test/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistanceTest.java @@ -256,7 +256,8 @@ private void checkHausdorff(String wkt1, String wkt2, String wktExpected) Geometry g1 = read(wkt1); Geometry g2 = read(wkt2); - Geometry result = DirectedHausdorffDistance.hausdorffDistanceLine(g1, g2); + Coordinate[] pts = DirectedHausdorffDistance.hausdorffDistancePoints(g1, g2); + Geometry result = g1.getFactory().createLineString(pts); Geometry expected = read(wktExpected); checkEqualExact(expected, result, TOLERANCE); } @@ -265,7 +266,8 @@ private void checkDistance(String wkt1, String wkt2, String wktExpected) { Geometry g1 = read(wkt1); Geometry g2 = read(wkt2); - Geometry result = DirectedHausdorffDistance.distanceLine(g1, g2); + Coordinate[] pts = DirectedHausdorffDistance.distancePoints(g1, g2); + Geometry result = g1.getFactory().createLineString(pts); Geometry expected = read(wktExpected); checkEqualExact(expected, result, TOLERANCE); } @@ -275,7 +277,8 @@ private void checkDistanceStartPtLen(String wkt1, String wkt2, Geometry g1 = read(wkt1); Geometry g2 = read(wkt2); - Geometry result = DirectedHausdorffDistance.distanceLine(g1, g2); + Coordinate[] pts = DirectedHausdorffDistance.distancePoints(g1, g2); + Geometry result = g1.getFactory().createLineString(pts); Geometry expected = read(wktExpected); Coordinate resultPt = result.getCoordinates()[0]; @@ -292,9 +295,7 @@ private void checkDistance(String wkt1, String wkt2, double tolerance, Geometry g1 = read(wkt1); Geometry g2 = read(wkt2); - Geometry result = DirectedHausdorffDistance.distanceLine(g1, g2, tolerance); - - double distResult = result.getLength(); + double distResult = DirectedHausdorffDistance.distance(g1, g2, tolerance); assertEquals(expectedDistance, distResult, TOLERANCE); } @@ -303,9 +304,7 @@ private void checkDistance(String wkt1, String wkt2, Geometry g1 = read(wkt1); Geometry g2 = read(wkt2); - Geometry result = DirectedHausdorffDistance.distanceLine(g1, g2); - - double distResult = result.getLength(); + double distResult = DirectedHausdorffDistance.distance(g1, g2); assertEquals(expectedDistance, distResult, TOLERANCE); } @@ -314,7 +313,6 @@ private void checkFullyWithinDistance(String a, String b, double distance, boole Geometry g2 = read(b); boolean result = DirectedHausdorffDistance.isFullyWithinDistance(g1, g2, distance); - assertEquals(expected, result); } diff --git a/modules/core/src/test/java/test/jts/perf/operation/distance/HausdorffDistancePerfTest.java b/modules/core/src/test/java/test/jts/perf/operation/distance/HausdorffDistancePerfTest.java index 9b4c113d5f..d1d0876d5d 100644 --- a/modules/core/src/test/java/test/jts/perf/operation/distance/HausdorffDistancePerfTest.java +++ b/modules/core/src/test/java/test/jts/perf/operation/distance/HausdorffDistancePerfTest.java @@ -14,10 +14,9 @@ import java.util.ArrayList; import java.util.List; -import org.locationtech.jts.algorithm.distance.DiscreteHausdorffDistance; import org.locationtech.jts.algorithm.distance.DirectedHausdorffDistance; +import org.locationtech.jts.algorithm.distance.DiscreteHausdorffDistance; import org.locationtech.jts.geom.Coordinate; -import org.locationtech.jts.geom.Envelope; import org.locationtech.jts.geom.Geometry; import org.locationtech.jts.geom.GeometryFactory; import org.locationtech.jts.geom.Point; @@ -71,9 +70,9 @@ public void runFullyWithin() { iter++; Debug.println(iter + " ----------------"); Debug.println(g); - //checkHausdorff(g, b); + checkHausdorff(g, b); //boolean isWithin = DirectedHausdorffDistance.isFullyWithinDistance(g, b, 4 * distance, 0.1); - boolean isWithin = dhd.isFullyWithinDistance(g, 4 * distance, 0.1); + boolean isWithin = dhd.isFullyWithinDistance(g, 4 * distance); //if (iter > 10) break; } @@ -82,7 +81,7 @@ public void runFullyWithin() { private void checkHausdorff(Geometry g, Geometry b) { double distDHD = DiscreteHausdorffDistance.distance(g, b, 0.01); - double distHD = DirectedHausdorffDistance.hausdorffDistanceLine(g, b).getLength(); + double distHD = DirectedHausdorffDistance.hausdorffDistance(g, b); //-- performance testing only //double distDHD = distHD; @@ -97,7 +96,7 @@ private void checkHausdorff(Geometry g, Geometry b) { } private void checkFullyWithinDistance(Geometry g, Geometry b) { - double distDHD = DirectedHausdorffDistance.distanceLine(g, b, 0.01).getLength(); + double distDHD = DirectedHausdorffDistance.distance(g, b, 0.01); double tol = distDHD / 1000; boolean isWithin = DirectedHausdorffDistance.isFullyWithinDistance(g, b, 1.05 * distDHD, tol); boolean isBeyond = ! DirectedHausdorffDistance.isFullyWithinDistance(g, b, 0.95 * distDHD, tol); diff --git a/modules/core/src/test/java/test/jts/perf/operation/distance/HausdorffDistanceStressTest.java b/modules/core/src/test/java/test/jts/perf/operation/distance/HausdorffDistanceStressTest.java index 9fa4d1b98e..3b7fb1693b 100644 --- a/modules/core/src/test/java/test/jts/perf/operation/distance/HausdorffDistanceStressTest.java +++ b/modules/core/src/test/java/test/jts/perf/operation/distance/HausdorffDistanceStressTest.java @@ -48,6 +48,7 @@ private void run() { Debug.println(iter + " ----------------"); Debug.println(g); checkHausdorff(g, b); + checkFullyWithinDistance(g, b); //if (iter > 10) break; } @@ -56,7 +57,7 @@ private void run() { private void checkHausdorff(Geometry g, Geometry b) { double distDHD = DiscreteHausdorffDistance.distance(g, b, 0.01); - double distHD = DirectedHausdorffDistance.hausdorffDistanceLine(g, b).getLength(); + double distHD = DirectedHausdorffDistance.hausdorffDistance(g, b); //-- performance testing only //double distDHD = distHD; @@ -66,19 +67,18 @@ private void checkHausdorff(Geometry g, Geometry b) { if (err > .01) { System.out.println("<<<<<<<<<<< ERROR!"); } - - checkFullyWithinDistance(g, b); } private void checkFullyWithinDistance(Geometry g, Geometry b) { - double distDHD = DirectedHausdorffDistance.distanceLine(g, b, 0.01).getLength(); - double tol = distDHD / 1000; - boolean isWithin = DirectedHausdorffDistance.isFullyWithinDistance(g, b, 1.05 * distDHD, tol); - boolean isBeyond = ! DirectedHausdorffDistance.isFullyWithinDistance(g, b, 0.95 * distDHD, tol); + double distDHD = DirectedHausdorffDistance.distance(g, b, 0.01); + + boolean isWithin = DirectedHausdorffDistance.isFullyWithinDistance(g, b, 1.05 * distDHD); + boolean isBeyond = ! DirectedHausdorffDistance.isFullyWithinDistance(g, b, 0.95 * distDHD); + if (! (isWithin && isBeyond)) { System.out.format("ioWithin = %b isBeyond = %b\n", isWithin, isBeyond); System.out.println("<<<<<<<<<<< ERROR!"); - DirectedHausdorffDistance.isFullyWithinDistance(g, b, 0.75 * distDHD,tol); + //DirectedHausdorffDistance.isFullyWithinDistance(g, b, 0.75 * distDHD,tol); } } From d3235acbf30cba8074faa4733b540988487644fc Mon Sep 17 00:00:00 2001 From: Martin Davis Date: Mon, 23 Feb 2026 22:48:23 -0800 Subject: [PATCH 37/43] clean up Hausdorff distance functions --- .../jtstest/function/DistanceFunctions.java | 14 ++++---------- 1 file changed, 4 insertions(+), 10 deletions(-) diff --git a/modules/app/src/main/java/org/locationtech/jtstest/function/DistanceFunctions.java b/modules/app/src/main/java/org/locationtech/jtstest/function/DistanceFunctions.java index 10d7aa226b..e26b35335a 100644 --- a/modules/app/src/main/java/org/locationtech/jtstest/function/DistanceFunctions.java +++ b/modules/app/src/main/java/org/locationtech/jtstest/function/DistanceFunctions.java @@ -46,12 +46,6 @@ public static Geometry frechetDistanceLine(Geometry a, Geometry b) return a.getFactory().createLineString(dist.getCoordinates()); } - @Metadata(description="Oriented discrete Hausdorff distance line from A to B") - public static Geometry orientedDiscreteHausdorffLine(Geometry a, Geometry b) - { - return DiscreteHausdorffDistance.orientedDistanceLine(a, b); - } - @Metadata(description="Oriented discrete Hausdorff distance from A to B") public static double orientedDiscreteHausdorffDistance(Geometry a, Geometry b) { @@ -66,12 +60,12 @@ public static Geometry orientedDiscreteHausdorffLineDensify(Geometry a, Geometry return DiscreteHausdorffDistance.orientedDistanceLine(a, b, frac); } - @Metadata(description="Clipped Oriented Hausdorff distance from A to B") - public static Geometry clippedOrientedHausdorffDistanceLine(Geometry a, Geometry b) + @Metadata(description="Clipped directed Hausdorff distance from A to B") + public static Geometry clippedDirectedHausdorffLine(Geometry a, Geometry b) { - //TODO: would this be more efficient done as part of DiscreteHausdorffDistance? Geometry clippedLine = LinearReferencingFunctions.project(a, b); - return DiscreteHausdorffDistance.orientedDistanceLine(clippedLine, b); + Coordinate[] pts = DirectedHausdorffDistance.distancePoints(clippedLine, b); + return a.getFactory().createLineString(pts); } @Metadata(description="Directed Hausdorff distance from A to B, up to tolerance") From 66c749d66f1ef6494eb359b569131cbde210fb78 Mon Sep 17 00:00:00 2001 From: Martin Davis Date: Tue, 24 Feb 2026 10:58:21 -0800 Subject: [PATCH 38/43] Handle empty inputs --- .../distance/DirectedHausdorffDistance.java | 62 ++++++++++++------- .../DirectedHausdorffDistanceTest.java | 52 ++++++++++++++++ 2 files changed, 92 insertions(+), 22 deletions(-) diff --git a/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java b/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java index cca9f4a5a4..a6529e9b14 100644 --- a/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java +++ b/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java @@ -30,13 +30,12 @@ import org.locationtech.jts.operation.distance.IndexedFacetDistance; /** - * Computes the directed Hausdorff distance from one geometry to another, - * up to an approximation distance tolerance. + * Computes the directed Hausdorff distance from one geometry to another. * The directed Hausdorff distance is the maximum distance any point * on a query geometry A can be from a target geometry B. * Equivalently, every point in the query geometry is within that distance * of the target geometry. - * The class can compute a pair of points at which the DHD is obtained: + * The class can compute a pair of points at which the distance is attained: * [ farthest A point, nearest B point ]. *

*The directed Hausdorff distance (DHD) is defined as: @@ -51,7 +50,7 @@ * HD(A,B) = max(DHD(A,B), DHD(B,A)) * * This can be computed via the - * {@link #hausdorffDistanceLine(Geometry, Geometry, double)} function. + * {@link #hausdorffDistancePoints(Geometry, Geometry)} function. *

* Points, lines and polygons are supported as input. * If the query geometry is polygonal, @@ -60,22 +59,22 @@ *

* A common use case is to test whether a geometry A lies fully within a given * distance of another one B. - * {@link #isFullyWithinDistance(Geometry, double, double)} + * {@link #isFullyWithinDistance(Geometry, double)} * can be used to test this efficiently. * It implements heuristic checks and short-circuiting to improve performance. - * This can much more efficient than computing whether A is covered by B.buffer(distance). - * It is also more accurate, since constructed buffers - * are only linearized approximations to the true buffer. *

* The class can be used in prepared mode. * Creating an instance on a target geometry caches indexes for that geometry. - * Then {@link #farthestPoints(Geometry, double) - * or {@link #isFullyWithinDistance(Geometry, double, double)} - * can be called efficiently for multiple query geometries. + * Then {@link #farthestPoints(Geometry) + * or {@link #isFullyWithinDistance(Geometry, double)} + * can be called efficiently for multiple query geometries. *

- * Due to the nature of the Hausdorff distance, - * performance is not very sensitive to the distance tolerance, - * so using a small tolerance is recommended. + * If the Hausdorff distance is attained at a non-vertex of the query geometry, + * the location is approximated. + * The algorithm uses a distance tolerance to control the approximation accuracy. + * The tolerance is automatically determined to balance between accuracy and performance. + * if more accuracy is desired some function signatures are provided + * which allow specifying a distance tolerance . *

* This algorithm is easier to use, more accurate, * and much faster than {@link DiscreteHausdorffDistance}. @@ -85,13 +84,16 @@ */ public class DirectedHausdorffDistance { + private static final double EMPTY_DISTANCE = Double.NaN; + /** * Computes the directed Hausdorff distance * of a query geometry A from a target one B. * * @param a the query geometry * @param b the target geometry - * @return the directed Hausdorff distance + * @return the directed Hausdorff distance, + * or NaN if an input is empty */ public static double distance(Geometry a, Geometry b) { @@ -107,7 +109,8 @@ public static double distance(Geometry a, Geometry b) * @param a the query geometry * @param b the target geometry * @param tolerance the accuracy distance tolerance - * @return the directed Hausdorff distance + * @return the directed Hausdorff distance, + * or NaN if an input is empty */ public static double distance(Geometry a, Geometry b, double tolerance) { @@ -122,7 +125,8 @@ public static double distance(Geometry a, Geometry b, double tolerance) * @param a the query geometry * @param b the target geometry * @param tolerance the accuracy distance tolerance - * @return a pair of points [ptA, ptB] demonstrating the distance + * @return a pair of points [ptA, ptB] demonstrating the distance, + * or null if an input is empty */ public static Coordinate[] distancePoints(Geometry a, Geometry b) { @@ -137,7 +141,8 @@ public static Coordinate[] distancePoints(Geometry a, Geometry b) * @param a the query geometry * @param b the target geometry * @param tolerance the accuracy distance tolerance - * @return a pair of points [ptA, ptB] demonstrating the distance + * @return a pair of points [ptA, ptB] demonstrating the distance, + * or null if an input is empty */ public static Coordinate[] distancePoints(Geometry a, Geometry b, double tolerance) { @@ -152,7 +157,8 @@ public static Coordinate[] distancePoints(Geometry a, Geometry b, double toleran * * @param a a geometry * @param b a geometry - * @return a pair of points [ptA, ptB] demonstrating the Hausdorff distance + * @return a pair of points [ptA, ptB] demonstrating the Hausdorff distance, + * or null if an input is empty */ public static Coordinate[] hausdorffDistancePoints(Geometry a, Geometry b) { @@ -176,7 +182,7 @@ public static Coordinate[] hausdorffDistancePoints(Geometry a, Geometry b) * * @param a a geometry * @param b a geometry - * @return the Hausdorff distance + * @return the Hausdorff distance, or NaN if an input is empty */ public static double hausdorffDistance(Geometry a, Geometry b) { @@ -220,6 +226,8 @@ public static boolean isFullyWithinDistance(Geometry a, Geometry b, double maxDi } private static double distance(Coordinate[] pts) { + if (pts == null) + return EMPTY_DISTANCE; return pts[0].distance(pts[1]); } @@ -289,6 +297,9 @@ public boolean isFullyWithinDistance(Geometry a, double maxDistance, double tole return false; Coordinate[] maxDistCoords = computeDistancePoints(a, tolerance, maxDistance); + //-- handle empty case + if (maxDistCoords == null) + return false; return distance(maxDistCoords) <= maxDistance; } @@ -316,9 +327,11 @@ private static boolean isBeyond(Envelope envA, Envelope envB, double maxDistance /** * Computes a pair of points which attain the directed Hausdorff distance * of a query geometry A from the target B. + * If either geometry is empty the result is null. * * @param geomA the query geometry - * @return a pair of points [ptA, ptB] attaining the distance + * @return a pair of points [ptA, ptB] attaining the distance, + * or null if an input is empty */ public Coordinate[] farthestPoints(Geometry geomA) { double tolerance = computeTolerance(geomA); @@ -329,16 +342,21 @@ public Coordinate[] farthestPoints(Geometry geomA) { * Computes a pair of points which attain the directed Hausdorff distance * of a query geometry A from the target B, * up to a given distance accuracy. + * If either geometry is empty the result is null. * * @param geomA the query geometry * @param tolerance the approximation distance tolerance - * @return a pair of points [ptA, ptB] attaining the distance + * @return a pair of points [ptA, ptB] attaining the distance, + * or null if an input is empty */ public Coordinate[] farthestPoints(Geometry geomA, double tolerance) { return computeDistancePoints(geomA, tolerance, -1.0); } private Coordinate[] computeDistancePoints(Geometry geomA, double tolerance, double maxDistanceLimit) { + if (geomA.isEmpty() || geomB.isEmpty()) + return null; + if (geomA.getDimension() == Dimension.P) { return computeForPoints(geomA, maxDistanceLimit); } diff --git a/modules/core/src/test/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistanceTest.java b/modules/core/src/test/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistanceTest.java index fa77d17670..175eca0cde 100644 --- a/modules/core/src/test/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistanceTest.java +++ b/modules/core/src/test/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistanceTest.java @@ -26,6 +26,22 @@ public static void main(String args[]) { public DirectedHausdorffDistanceTest(String name) { super(name); } + public void testEmptyPoint() + { + checkDistanceEmpty("POINT EMPTY", "POINT (1 1)"); + } + + public void testEmptyLine() + { + checkDistanceEmpty("LINESTRING EMPTY", "LINESTRING (0 0, 2 1)"); + } + + public void testEmptyPolygon() + { + checkDistanceEmpty("POLYGON EMPTY", "POLYGON ((1 9, 9 9, 9 1, 1 1, 1 9))"); + } + + //-------------------------------------- public void testPointPoint() { checkHausdorff("POINT (0 0)", "POINT (1 1)", @@ -191,6 +207,20 @@ public void testInteriorSegmentsSameExterior() //----------------------------------------------------- + public void testFullyWithinDistanceEmptyPoints() + { + String a = "POINT EMPTY"; + String b = "MULTIPOINT ((1 1), (9 9))"; + checkFullyWithinDistanceEmpty(a, b); + } + + public void testFullyWithinDistanceEmptyLine() + { + String a = "LINESTRING EMPTY"; + String b = "LINESTRING (9 9, 1 1)"; + checkFullyWithinDistanceEmpty(a, b); + } + //-- shows withinDistance envelope check not triggering for disconnected A public void testFullyWithinDistancePoints() { @@ -316,4 +346,26 @@ private void checkFullyWithinDistance(String a, String b, double distance, boole assertEquals(expected, result); } + private void checkFullyWithinDistanceEmpty(String a, String b) { + checkFullyWithinDistance(a, b, 0, false); + checkFullyWithinDistance(b, a, 0, false); + checkFullyWithinDistance(a, b, 1, false); + checkFullyWithinDistance(b, a, 1, false); + checkFullyWithinDistance(a, b, 1000, false); + checkFullyWithinDistance(b, a, 1000, false); + } + + private void checkDistanceEmpty(String a, String b) { + Geometry g1 = read(a); + Geometry g2 = read(b); + + Coordinate[] resultPts = DirectedHausdorffDistance.distancePoints(g1, g2); + assert(resultPts == null); + + double resultDist = DirectedHausdorffDistance.distance(g1, g2); + assert(Double.isNaN(resultDist)); + + double hausdorffDist = DirectedHausdorffDistance.hausdorffDistance(g1, g2); + assert(Double.isNaN(hausdorffDist)); + } } From 08b81e411c830757d2d4046b2fac7579f4ae9d38 Mon Sep 17 00:00:00 2001 From: Martin Davis Date: Tue, 24 Feb 2026 11:07:25 -0800 Subject: [PATCH 39/43] Variable renames --- .../distance/DirectedHausdorffDistance.java | 118 +++++++++--------- 1 file changed, 59 insertions(+), 59 deletions(-) diff --git a/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java b/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java index a6529e9b14..e313008aa5 100644 --- a/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java +++ b/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java @@ -249,17 +249,17 @@ private static double computeTolerance(Geometry geom) { return geom.getEnvelopeInternal().getDiameter() / AUTO_TOLERANCE_FACTOR; } - private Geometry geomB; - private TargetDistance distanceToB; + private Geometry target; + private TargetDistance targetDistance; /** - * Create a new instance for a target geometry B + * Create a new instance for a target geometry. * - * @param b the geometry to compute the distance from + * @param geom the geometry to compute the distance from */ - public DirectedHausdorffDistance(Geometry b) { - geomB = b; - distanceToB = new TargetDistance(geomB); + public DirectedHausdorffDistance(Geometry geom) { + this.target = geom; + targetDistance = new TargetDistance(this.target); } /** @@ -268,15 +268,15 @@ public DirectedHausdorffDistance(Geometry b) { * from the target than the specified distance. * This is the case if DHD(A, B) > maxDistance. * - * @param a the query geometry + * @param geom the query geometry * @param maxDistance the distance limit * @return true if the query geometry lies fully within the distance of the target */ - public boolean isFullyWithinDistance(Geometry a, double maxDistance) { + public boolean isFullyWithinDistance(Geometry geom, double maxDistance) { //TODO: should the tolerance be computed as a fraction of the maxDistance? //return isFullyWithinDistance(a, maxDistance, computeTolerance(a)); double tolerance = maxDistance / AUTO_TOLERANCE_FACTOR; - return isFullyWithinDistance(a, maxDistance, tolerance); + return isFullyWithinDistance(geom, maxDistance, tolerance); } /** @@ -286,17 +286,17 @@ public boolean isFullyWithinDistance(Geometry a, double maxDistance) { * from the target than the specified distance. * This is the case if DHD(A, B) > maxDistance. * - * @param a the query geometry + * @param geom the query geometry * @param maxDistance the distance limit * @param tolerance the accuracy distance tolerance * @return true if the query geometry lies fully within the distance of the target */ - public boolean isFullyWithinDistance(Geometry a, double maxDistance, double tolerance) { + public boolean isFullyWithinDistance(Geometry geom, double maxDistance, double tolerance) { //-- envelope checks - if (isBeyond(a.getEnvelopeInternal(), geomB.getEnvelopeInternal(), maxDistance)) + if (isBeyond(geom.getEnvelopeInternal(), target.getEnvelopeInternal(), maxDistance)) return false; - Coordinate[] maxDistCoords = computeDistancePoints(a, tolerance, maxDistance); + Coordinate[] maxDistCoords = computeDistancePoints(geom, tolerance, maxDistance); //-- handle empty case if (maxDistCoords == null) return false; @@ -304,13 +304,13 @@ public boolean isFullyWithinDistance(Geometry a, double maxDistance, double tole } /** - * Tests if a geometry must have a point farther than the maximum distance - * using the geometry envelopes. + * Tests if an envelope must have a side farther than the maximum distance + * from another envelope. * * @param envA * @param envB * @param maxDistance - * @return true if geometry A must have a far point from B + * @return true if envelope A must have a far point from envelope B */ private static boolean isBeyond(Envelope envA, Envelope envB, double maxDistance) { /** @@ -329,13 +329,13 @@ private static boolean isBeyond(Envelope envA, Envelope envB, double maxDistance * of a query geometry A from the target B. * If either geometry is empty the result is null. * - * @param geomA the query geometry + * @param geom the query geometry * @return a pair of points [ptA, ptB] attaining the distance, * or null if an input is empty */ - public Coordinate[] farthestPoints(Geometry geomA) { - double tolerance = computeTolerance(geomA); - return farthestPoints(geomA, tolerance); + public Coordinate[] farthestPoints(Geometry geom) { + double tolerance = computeTolerance(geom); + return farthestPoints(geom, tolerance); } /** @@ -344,24 +344,24 @@ public Coordinate[] farthestPoints(Geometry geomA) { * up to a given distance accuracy. * If either geometry is empty the result is null. * - * @param geomA the query geometry + * @param geom the query geometry * @param tolerance the approximation distance tolerance * @return a pair of points [ptA, ptB] attaining the distance, * or null if an input is empty */ - public Coordinate[] farthestPoints(Geometry geomA, double tolerance) { - return computeDistancePoints(geomA, tolerance, -1.0); + public Coordinate[] farthestPoints(Geometry geom, double tolerance) { + return computeDistancePoints(geom, tolerance, -1.0); } - private Coordinate[] computeDistancePoints(Geometry geomA, double tolerance, double maxDistanceLimit) { - if (geomA.isEmpty() || geomB.isEmpty()) + private Coordinate[] computeDistancePoints(Geometry geom, double tolerance, double maxDistanceLimit) { + if (geom.isEmpty() || target.isEmpty()) return null; - if (geomA.getDimension() == Dimension.P) { - return computeForPoints(geomA, maxDistanceLimit); + if (geom.getDimension() == Dimension.P) { + return computeForPoints(geom, maxDistanceLimit); } //TODO: handle mixed geoms with points - Coordinate[] maxDistPtsEdge = computeForEdges(geomA, tolerance, maxDistanceLimit); + Coordinate[] maxDistPtsEdge = computeForEdges(geom, tolerance, maxDistanceLimit); if (isBeyondLimit(distance(maxDistPtsEdge), maxDistanceLimit)) { return maxDistPtsEdge; @@ -370,8 +370,8 @@ private Coordinate[] computeDistancePoints(Geometry geomA, double tolerance, dou /** * Polygonal query geometry may have an interior point as the farthest point. */ - if (geomA.getDimension() == Dimension.A) { - Coordinate[] maxDistPtsInterior = computeForAreaInterior(geomA, tolerance); + if (geom.getDimension() == Dimension.A) { + Coordinate[] maxDistPtsInterior = computeForAreaInterior(geom, tolerance); if (maxDistPtsInterior != null && distance(maxDistPtsInterior) > distance(maxDistPtsEdge)) { return maxDistPtsInterior; @@ -380,20 +380,20 @@ && distance(maxDistPtsInterior) > distance(maxDistPtsEdge)) { return maxDistPtsEdge; } - private Coordinate[] computeForPoints(Geometry geomA, double maxDistanceLimit) { + private Coordinate[] computeForPoints(Geometry geom, double maxDistanceLimit) { double maxDist = -1.0;; Coordinate[] maxDistPtsAB = null; - Iterator geomi = new GeometryCollectionIterator(geomA); + Iterator geomi = new GeometryCollectionIterator(geom); while (geomi.hasNext()) { - Geometry geomElemA = (Geometry) geomi.next(); - if (! (geomElemA instanceof Point)) + Geometry geomElem = (Geometry) geomi.next(); + if (! (geomElem instanceof Point)) continue; - Coordinate pA = geomElemA.getCoordinate(); - Coordinate pB = distanceToB.nearestPoint(pA); + Coordinate pA = geomElem.getCoordinate(); + Coordinate pB = targetDistance.nearestPoint(pA); double dist = pA.distance(pB); - boolean isInterior = dist > 0 && distanceToB.isInterior(pA); + boolean isInterior = dist > 0 && targetDistance.isInterior(pA); //-- check for interior point if (isInterior) { dist = 0; @@ -410,8 +410,8 @@ private Coordinate[] computeForPoints(Geometry geomA, double maxDistanceLimit) { return maxDistPtsAB; } - private Coordinate[] computeForEdges(Geometry geomA, double tolerance, double maxDistanceLimit) { - PriorityQueue segQueue = createSegQueue(geomA); + private Coordinate[] computeForEdges(Geometry geom, double tolerance, double maxDistanceLimit) { + PriorityQueue segQueue = createSegQueue(geom); DHDSegment segMaxDist = null; long iter = 0; @@ -476,7 +476,7 @@ private Coordinate[] computeForEdges(Geometry geomA, double tolerance, double ma * so bisect and keep searching */ if ((segMaxBound.getLength() > tolerance)) { - DHDSegment[] bisects = segMaxBound.bisect(distanceToB); + DHDSegment[] bisects = segMaxBound.bisect(targetDistance); addNonInterior(bisects[0], segQueue); addNonInterior(bisects[1], segQueue); } @@ -494,13 +494,13 @@ private Coordinate[] computeForEdges(Geometry geomA, double tolerance, double ma * In this case distance is zero. * Return a single coordinate of the input as a representative point */ - Coordinate maxPt = geomA.getCoordinate(); + Coordinate maxPt = geom.getCoordinate(); return pair(maxPt, maxPt); } private boolean isSameOrCollinear(DHDSegment seg) { - CoordinateSequenceLocation f0 = distanceToB.nearestLocation(seg.p0); - CoordinateSequenceLocation f1 = distanceToB.nearestLocation(seg.p1); + CoordinateSequenceLocation f0 = targetDistance.nearestLocation(seg.p0); + CoordinateSequenceLocation f1 = targetDistance.nearestLocation(seg.p1); return f0.isSameSegment(f1); } @@ -530,11 +530,11 @@ private boolean isInterior(DHDSegment segment) { if (segment.getMaxDistance() > 0.0) { return false; } - return distanceToB.isInterior(segment.getEndpoint(0), segment.getEndpoint(1)); + return targetDistance.isInterior(segment.getEndpoint(0), segment.getEndpoint(1)); } /** - * If the query geometry A is polygonal, it is possible + * If the query geometry is polygonal, it is possible * the farthest point lies in its interior. * In this case it occurs at the centre of the Largest Empty Circle * with B as obstacles and query geometry as constraint. @@ -545,20 +545,20 @@ private boolean isInterior(DHDSegment segment) { * can determine that the max distance is at * the previously computed edge points. * - * @param geomA + * @param geom * @param tolerance - * @return the maximum distance point pair at an interior point of A, + * @return the maximum distance point pair at an interior point of a polygon, * or null if it is known to not occur at an interior point */ - private Coordinate[] computeForAreaInterior(Geometry geomA, double tolerance) { + private Coordinate[] computeForAreaInterior(Geometry geom, double tolerance) { //TODO: extract polygonal geoms from A - Geometry polygonalA = geomA; + Geometry polygonal = geom; /** * Optimization - skip if A interior cannot intersect B, * and thus farther point must lie on A boundary */ - if (polygonalA.getEnvelopeInternal().disjoint(geomB.getEnvelopeInternal())) { + if (polygonal.getEnvelopeInternal().disjoint(target.getEnvelopeInternal())) { return null; } @@ -572,17 +572,17 @@ private Coordinate[] computeForAreaInterior(Geometry geomA, double tolerance) { //TODO: add short-circuiting based on maxDistanceLimit? - Point centerPt = LargestEmptyCircle.getCenter(geomB, polygonalA, lecTol); + Point centerPt = LargestEmptyCircle.getCenter(target, polygonal, lecTol); Coordinate ptA = centerPt.getCoordinate(); /** * If LEC centre is in B, the max distance is zero, so return null. * This will cause the computed segment distance to be returned, * which is preferred since it occurs on a vertex or edge. */ - if (distanceToB.isInterior(ptA)) { + if (targetDistance.isInterior(ptA)) { return null; } - Coordinate ptB = distanceToB.nearestFacetPoint(ptA); + Coordinate ptB = targetDistance.nearestFacetPoint(ptA); return pair(ptA, ptB); } @@ -593,12 +593,12 @@ private Coordinate[] computeForAreaInterior(Geometry geomA, double tolerance) { * So if a query geometry is fully covered by the target * the returned queue is empty. * - * @param geomA + * @param geom * @return the segment priority queue */ - private PriorityQueue createSegQueue(Geometry geomA) { + private PriorityQueue createSegQueue(Geometry geom) { PriorityQueue priq = new PriorityQueue(); - geomA.apply(new GeometryComponentFilter() { + geom.apply(new GeometryComponentFilter() { @Override public void filter(Geometry geom) { @@ -622,11 +622,11 @@ private void addSegments(Coordinate[] pts, PriorityQueue priq) { for (int i = 0; i < pts.length - 1; i++) { DHDSegment seg; if (i == 0) { - seg = DHDSegment.create(pts[i], pts[i + 1], distanceToB); + seg = DHDSegment.create(pts[i], pts[i + 1], targetDistance); } else { //-- avoiding recomputing prev pt distance - seg = DHDSegment.create(prevSeg, pts[i + 1], distanceToB); + seg = DHDSegment.create(prevSeg, pts[i + 1], targetDistance); } prevSeg = seg; From fe5471cf0bf4df6ff57ba99716d0bff73213fa39 Mon Sep 17 00:00:00 2001 From: Martin Davis Date: Tue, 24 Feb 2026 13:49:19 -0800 Subject: [PATCH 40/43] Javadoc --- .../algorithm/distance/DirectedHausdorffDistance.java | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java b/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java index e313008aa5..51ba240ee3 100644 --- a/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java +++ b/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java @@ -57,10 +57,10 @@ * the point at maximum distance may occur in the interior of a query polygon. * For a polygonal target geometry the point always lies on the boundary. *

- * A common use case is to test whether a geometry A lies fully within a given - * distance of another one B. - * {@link #isFullyWithinDistance(Geometry, double)} - * can be used to test this efficiently. + * The directed Hausdorff distance can be used to test + * whether a geometry lies fully within a given distance of another one. + * The {@link #isFullyWithinDistance(Geometry, double)} function + * is provided to execute this test efficiently. * It implements heuristic checks and short-circuiting to improve performance. *

* The class can be used in prepared mode. @@ -70,7 +70,7 @@ * can be called efficiently for multiple query geometries. *

* If the Hausdorff distance is attained at a non-vertex of the query geometry, - * the location is approximated. + * the location must be approximated. * The algorithm uses a distance tolerance to control the approximation accuracy. * The tolerance is automatically determined to balance between accuracy and performance. * if more accuracy is desired some function signatures are provided From 6e2f202e068bdfe302a68f34e280c4235f84cd52 Mon Sep 17 00:00:00 2001 From: Martin Davis Date: Tue, 24 Feb 2026 14:22:45 -0800 Subject: [PATCH 41/43] Update debugging code --- .../algorithm/distance/DirectedHausdorffDistance.java | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java b/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java index 51ba240ee3..5350555f8b 100644 --- a/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java +++ b/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java @@ -423,20 +423,21 @@ private Coordinate[] computeForEdges(Geometry geom, double tolerance, double max //System.out.println(WKTWriter.toLineString(segMaxBound.getMaxDistPts())); /* - double maxDistBound = segMaxBound.maxDistanceBound(); - double maxDist = segMaxBound.maxDistance(); + double maxDistBound = segMaxBound.getMaxDistanceBound(); + double maxDist = segMaxBound.getMaxDistance(); System.out.format("%s len: %f bound: %f maxDist: %f\n", - segMaxBound, segMaxBound.length(), maxDistBound, maxDist); + segMaxBound, segMaxBound.getLength(), maxDistBound, maxDist); if (maxDist > 1) { System.out.println("FOUND"); } -*/ +//*/ /** * Save if segment point is farther than current farthest */ if (segMaxDist == null || segMaxBound.getMaxDistance() > segMaxDist.getMaxDistance()) { segMaxDist = segMaxBound; + //System.out.println(">>>> " + WKTWriter.toLineString(segMaxDist.getMaxDistPts())); } /** * Stop searching if remaining items in queue must all be closer From 63d79af380780ce057e01a97f33ce82271dbce1f Mon Sep 17 00:00:00 2001 From: Martin Davis Date: Tue, 24 Feb 2026 20:21:44 -0800 Subject: [PATCH 42/43] Refactor tolerance factors --- .../distance/DirectedHausdorffDistance.java | 29 ++++++++++++------- 1 file changed, 18 insertions(+), 11 deletions(-) diff --git a/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java b/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java index 5350555f8b..e5f9562813 100644 --- a/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java +++ b/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java @@ -83,8 +83,6 @@ * */ public class DirectedHausdorffDistance { - - private static final double EMPTY_DISTANCE = Double.NaN; /** * Computes the directed Hausdorff distance @@ -235,15 +233,28 @@ private static Coordinate[] pair(Coordinate p0, Coordinate p1) { return new Coordinate[] { p0.copy(), p1.copy() }; } + private static final double EMPTY_DISTANCE = Double.NaN; + + /** + * Heuristic automatic tolerance factor + */ + private static final double AUTO_TOLERANCE_FACTOR = 1.0e4; + /** * Heuristic factor to improve performance of area interior farthest point computation. + * The LargestEmptyCircle computation is much slower than the boundary one, + * is unlikely to occur, and accuracy is less critical (and obvious). + * Improve performance by using a coarser distance tolerance. */ - private static final double AREA_INTERIOR_PERFORMANCE_FACTOR = 20; + private static final double AREA_INTERIOR_TOLERANCE_FACTOR = 20; /** - * Heuristic automatic tolerance factor + * Tolerance factor for isFullyWithinDistance. + * A larger factor is used to increase accuracy. + * The operation will usually short-circuit, so performance impact is low. */ - private static final double AUTO_TOLERANCE_FACTOR = 1.0e4; + private static final double FULLY_WITHIN_TOLERANCE_FACTOR = 10 * AUTO_TOLERANCE_FACTOR; + private static double computeTolerance(Geometry geom) { return geom.getEnvelopeInternal().getDiameter() / AUTO_TOLERANCE_FACTOR; @@ -273,9 +284,7 @@ public DirectedHausdorffDistance(Geometry geom) { * @return true if the query geometry lies fully within the distance of the target */ public boolean isFullyWithinDistance(Geometry geom, double maxDistance) { - //TODO: should the tolerance be computed as a fraction of the maxDistance? - //return isFullyWithinDistance(a, maxDistance, computeTolerance(a)); - double tolerance = maxDistance / AUTO_TOLERANCE_FACTOR; + double tolerance = maxDistance / FULLY_WITHIN_TOLERANCE_FACTOR; return isFullyWithinDistance(geom, maxDistance, tolerance); } @@ -569,11 +578,9 @@ private Coordinate[] computeForAreaInterior(Geometry geom, double tolerance) { * and accuracy is probably less critical (or obvious). * So improve performance by using a coarser distance tolerance. */ - double lecTol = AREA_INTERIOR_PERFORMANCE_FACTOR * tolerance; - //TODO: add short-circuiting based on maxDistanceLimit? - Point centerPt = LargestEmptyCircle.getCenter(target, polygonal, lecTol); + Point centerPt = LargestEmptyCircle.getCenter(target, polygonal, tolerance * AREA_INTERIOR_TOLERANCE_FACTOR); Coordinate ptA = centerPt.getCoordinate(); /** * If LEC centre is in B, the max distance is zero, so return null. From 50c0104a8e9e00bb80c9f15f23463bdd05b5eec6 Mon Sep 17 00:00:00 2001 From: Martin Davis Date: Fri, 27 Feb 2026 10:56:38 -0800 Subject: [PATCH 43/43] Add checks for invalid and zero input --- .../distance/DirectedHausdorffDistance.java | 38 +++++++-- .../DirectedHausdorffDistanceTest.java | 82 +++++++++++++++++++ 2 files changed, 114 insertions(+), 6 deletions(-) diff --git a/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java b/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java index e5f9562813..285a077fa6 100644 --- a/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java +++ b/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java @@ -255,7 +255,14 @@ private static Coordinate[] pair(Coordinate p0, Coordinate p1) { */ private static final double FULLY_WITHIN_TOLERANCE_FACTOR = 10 * AUTO_TOLERANCE_FACTOR; - + /** + * Computes an automatic distance tolerance for a query geometry. + * A tolerance of zero may be returned (e.g. for a zero-length line); + * this is supported by the distance algorithm. + * + * @param geom the query geometry + * @return the distance tolerance + */ private static double computeTolerance(Geometry geom) { return geom.getEnvelopeInternal().getDiameter() / AUTO_TOLERANCE_FACTOR; } @@ -363,6 +370,14 @@ public Coordinate[] farthestPoints(Geometry geom, double tolerance) { } private Coordinate[] computeDistancePoints(Geometry geom, double tolerance, double maxDistanceLimit) { + /** + * Negative tolerances are not allowed. + * Zero tolerance is allowed, to support zero-size input. + */ + if (tolerance < 0.0 ) { + throw new IllegalArgumentException("Tolerance must be non-negative"); + } + if (geom.isEmpty() || target.isEmpty()) return null; @@ -412,7 +427,8 @@ private Coordinate[] computeForPoints(Geometry geom, double maxDistanceLimit) { maxDist = dist; maxDistPtsAB = pair(pA, pB); } - if (isBeyondLimit(maxDist, maxDistanceLimit)) { + if (isValidLimit(maxDistanceLimit) + && isBeyondLimit(maxDist, maxDistanceLimit)) { break; } } @@ -460,10 +476,12 @@ private Coordinate[] computeForEdges(Geometry geom, double tolerance, double max * - if segment distance bound is less than distance limit, no other segment can be farther * - if a point of segment is farther than limit, isFulyWithin must be false */ - if (isWithinLimit(segMaxBound.getMaxDistanceBound(), maxDistanceLimit) + if (isValidLimit(maxDistanceLimit)) { + if (isWithinLimit(segMaxBound.getMaxDistanceBound(), maxDistanceLimit) || isBeyondLimit(segMaxBound.getMaxDistance(), maxDistanceLimit) ) { - break; + break; + } } /** @@ -483,9 +501,10 @@ private Coordinate[] computeForEdges(Geometry geom, double tolerance, double max /** * If segment is longer than tolerance * it might provide a better max distance point, - * so bisect and keep searching + * so bisect and keep searching. */ - if ((segMaxBound.getLength() > tolerance)) { + if (tolerance > 0 + && segMaxBound.getLength() > tolerance) { DHDSegment[] bisects = segMaxBound.bisect(targetDistance); addNonInterior(bisects[0], segQueue); addNonInterior(bisects[1], segQueue); @@ -514,6 +533,10 @@ private boolean isSameOrCollinear(DHDSegment seg) { return f0.isSameSegment(f1); } + private static boolean isValidLimit(double limit) { + return limit >= 0.0; + } + private static boolean isBeyondLimit(double maxDist, double maxDistanceLimit) { return maxDistanceLimit >= 0 && maxDist > maxDistanceLimit; } @@ -561,6 +584,9 @@ private boolean isInterior(DHDSegment segment) { * or null if it is known to not occur at an interior point */ private Coordinate[] computeForAreaInterior(Geometry geom, double tolerance) { + if (tolerance <= 0.0) + return null; + //TODO: extract polygonal geoms from A Geometry polygonal = geom; diff --git a/modules/core/src/test/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistanceTest.java b/modules/core/src/test/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistanceTest.java index 175eca0cde..2e83ef0ff3 100644 --- a/modules/core/src/test/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistanceTest.java +++ b/modules/core/src/test/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistanceTest.java @@ -26,6 +26,8 @@ public static void main(String args[]) { public DirectedHausdorffDistanceTest(String name) { super(name); } + //-- empty inputs + public void testEmptyPoint() { checkDistanceEmpty("POINT EMPTY", "POINT (1 1)"); @@ -42,6 +44,76 @@ public void testEmptyPolygon() } //-------------------------------------- + //-- extreme and invalid inputs + + public void testZeroTolerancePoint() + { + checkDistance("POINT (5 5)", "LINESTRING (5 1, 9 5)", + 0, + "LINESTRING (5 5, 7 3)"); + } + + public void testZeroToleranceLine() + { + checkDistance("LINESTRING (1 5, 5 5)", "LINESTRING (5 1, 9 5)", + 0, + "LINESTRING (1 5, 5 1)"); + } + + public void testZeroToleranceZeroLengthLineQuery() + { + checkDistance("LINESTRING (5 5, 5 5)", "LINESTRING (5 1, 9 5)", + 0, + "LINESTRING (5 5, 7 3)"); + } + + public void testZeroLengthLineQuery() + { + checkDistance("LINESTRING (5 5, 5 5)", "LINESTRING (5 1, 9 5)", + "LINESTRING (5 5, 7 3)"); + } + + public void testZeroLengthPolygonQuery() + { + checkDistance("POLYGON ((5 5, 5 5, 5 5, 5 5))", "LINESTRING (5 1, 9 5)", + "LINESTRING (5 5, 7 3)"); + } + + public void testZeroLengthLineTarget() + { + checkDistance("POINT (5 5)", "LINESTRING (5 1, 5 1)", + "LINESTRING (5 5, 5 1)"); + } + + public void testNegativeTolerancePoint() + { + try { + checkDistance("POINT (5 5)", "LINESTRING (5 1, 9 5)", + -1, + "LINESTRING (5 5, 7 3)"); + fail(); + } + catch (IllegalArgumentException expected) { + + } + } + + public void testNegativeToleranceLine() + { + try { + checkDistance("LINESTRING (1 5, 5 5)", "LINESTRING (5 1, 9 5)", + -1, + "LINESTRING (1 5, 5 1)"); + fail(); + } + catch (IllegalArgumentException expected) { + + } + } + + //-------------------------------------- + + public void testPointPoint() { checkHausdorff("POINT (0 0)", "POINT (1 1)", @@ -302,6 +374,16 @@ private void checkDistance(String wkt1, String wkt2, String wktExpected) { checkEqualExact(expected, result, TOLERANCE); } + private void checkDistance(String wkt1, String wkt2, double tolerance, String wktExpected) { + Geometry g1 = read(wkt1); + Geometry g2 = read(wkt2); + + Coordinate[] pts = DirectedHausdorffDistance.distancePoints(g1, g2, tolerance); + Geometry result = g1.getFactory().createLineString(pts); + Geometry expected = read(wktExpected); + checkEqualExact(expected, result, TOLERANCE); + } + private void checkDistanceStartPtLen(String wkt1, String wkt2, String wktExpected, double resultTolerance) { Geometry g1 = read(wkt1);