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..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 @@ -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; @@ -45,53 +46,61 @@ 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 from A to B") + public static double orientedDiscreteHausdorffDistance(Geometry a, Geometry b) + { + return DiscreteHausdorffDistance.orientedDistance(a, 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.distance(a, b); + return DiscreteHausdorffDistance.orientedDistanceLine(a, b, frac); + } + + @Metadata(description="Clipped directed Hausdorff distance from A to B") + public static Geometry clippedDirectedHausdorffLine(Geometry a, Geometry b) + { + Geometry clippedLine = LinearReferencingFunctions.project(a, b); + Coordinate[] pts = DirectedHausdorffDistance.distancePoints(clippedLine, b); + return a.getFactory().createLineString(pts); } - @Metadata(description="Hausdorff distance between A and B") - public static Geometry hausdorffDistanceLine(Geometry a, Geometry b) + @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.distanceLine(a, b); + return DirectedHausdorffDistance.distance(a, b, distTol); } - - @Metadata(description="Hausdorff distance between A and B, densified") - public static Geometry hausdorffDistanceLineDensify(Geometry a, Geometry b, - @Metadata(title="Densify fraction") - double frac) - { - return DiscreteHausdorffDistance.distanceLine(a, b, frac); - } - - @Metadata(description="Oriented Hausdorff distance from A to B") - public static Geometry orientedHausdorffDistanceLine(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 DiscreteHausdorffDistance.orientedDistanceLine(a, b); + Coordinate[] pts = DirectedHausdorffDistance.distancePoints(a, b, distTol); + return a.getFactory().createLineString(pts); } - - @Metadata(description="Oriented Hausdorff distance from A to B") - public static Geometry clippedOrientedHausdorffDistanceLine(Geometry a, Geometry b) + + @Metadata(description="Directed Hausdorff distance line from A to B") + public static Geometry directedHausdorffLine(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(a, b); + return a.getFactory().createLineString(pts); } - - @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="Hausdorff distance between A and B, up to tolerance") + public static Geometry hausdorffLine(Geometry a, Geometry b) { - return DiscreteHausdorffDistance.orientedDistanceLine(a, b, frac); + Coordinate[] pts = DirectedHausdorffDistance.hausdorffDistancePoints(a, b); + return a.getFactory().createLineString(pts); } - + + //-------------------------------------------- + public static double distanceIndexed(Geometry a, Geometry b) { return IndexedFacetDistance.distance(a, b); } @@ -117,4 +126,5 @@ public static Geometry nearestPointsIndexedEachB(Geometry a, Geometry b) { return a.getFactory().createMultiLineString(lines); } + } 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..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 @@ -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,25 @@ 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) { + return DirectedHausdorffDistance.isFullyWithinDistance(g, mask, maximumDistance); + } + }); + } + + 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) { + return dhd.isFullyWithinDistance(g, maximumDistance); + } + }); + } + 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/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 new file mode 100644 index 0000000000..285a077fa6 --- /dev/null +++ b/modules/core/src/main/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistance.java @@ -0,0 +1,844 @@ +/* + * 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.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.Envelope; +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.io.WKTWriter; +import org.locationtech.jts.operation.distance.CoordinateSequenceLocation; +import org.locationtech.jts.operation.distance.IndexedFacetDistance; + +/** + * 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 distance is attained: + * [ farthest A point, nearest B point ]. + *

+ *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 defined as: + * 

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

+ * Points, lines and polygons are supported as input. + * If the query geometry is polygonal, + * 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 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. + * Creating an instance on a target geometry caches indexes for that geometry. + * Then {@link #farthestPoints(Geometry) + * or {@link #isFullyWithinDistance(Geometry, double)} + * can be called efficiently for multiple query geometries. + *

+ * If the Hausdorff distance is attained at a non-vertex of the query geometry, + * 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 + * which allow specifying a distance tolerance . + *

+ * This algorithm is easier to use, more accurate, + * and much faster than {@link DiscreteHausdorffDistance}. + * + * @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 + * @return the directed Hausdorff distance, + * or NaN if an input is empty + */ + public static double distance(Geometry a, Geometry b) + { + DirectedHausdorffDistance hd = new DirectedHausdorffDistance(b); + return distance(hd.farthestPoints(a)); + } + + /** + * 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 accuracy distance tolerance + * @return the directed Hausdorff distance, + * or NaN if an input is empty + */ + public static double distance(Geometry a, Geometry b, double tolerance) + { + DirectedHausdorffDistance hd = new DirectedHausdorffDistance(b); + 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. + * + * @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, + * or null if an input is empty + */ + public static Coordinate[] distancePoints(Geometry a, Geometry b) + { + 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, 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, + * or null if an input is empty + */ + public static Coordinate[] distancePoints(Geometry a, Geometry b, double tolerance) + { + DirectedHausdorffDistance dhd = new DirectedHausdorffDistance(b); + return dhd.farthestPoints(a, 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 + * @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) + { + DirectedHausdorffDistance hdAB = new DirectedHausdorffDistance(b); + Coordinate[] ptsAB = hdAB.farthestPoints(a); + DirectedHausdorffDistance hdBA = new DirectedHausdorffDistance(a); + Coordinate[] ptsBA = hdBA.farthestPoints(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 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, or NaN if an input is empty + */ + public static double hausdorffDistance(Geometry a, Geometry b) + { + return distance(hausdorffDistancePoints(a, b)); + } + + /** + * 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 + * @return true if the query geometry lies fully within the distance of the target + */ + public static boolean isFullyWithinDistance(Geometry a, Geometry b, double maxDistance) + { + DirectedHausdorffDistance hd = new DirectedHausdorffDistance(b); + return hd.isFullyWithinDistance(a, 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, tolerance); + } + + private static double distance(Coordinate[] pts) { + if (pts == null) + return EMPTY_DISTANCE; + return pts[0].distance(pts[1]); + } + + 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_TOLERANCE_FACTOR = 20; + + /** + * 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 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; + } + + private Geometry target; + private TargetDistance targetDistance; + + /** + * Create a new instance for a target geometry. + * + * @param geom the geometry to compute the distance from + */ + public DirectedHausdorffDistance(Geometry geom) { + this.target = geom; + targetDistance = new TargetDistance(this.target); + } + + /** + * 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 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 geom, double maxDistance) { + double tolerance = maxDistance / FULLY_WITHIN_TOLERANCE_FACTOR; + return isFullyWithinDistance(geom, maxDistance, tolerance); + } + + /** + * 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 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 geom, double maxDistance, double tolerance) { + //-- envelope checks + if (isBeyond(geom.getEnvelopeInternal(), target.getEnvelopeInternal(), maxDistance)) + return false; + + Coordinate[] maxDistCoords = computeDistancePoints(geom, tolerance, maxDistance); + //-- handle empty case + if (maxDistCoords == null) + return false; + return distance(maxDistCoords) <= maxDistance; + } + + /** + * 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 envelope A must have a far point from envelope 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 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 + || envA.getMaxX() > envB.getMaxX() + maxDistance + || envA.getMaxY() > envB.getMaxY() + 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 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 geom) { + double tolerance = computeTolerance(geom); + return farthestPoints(geom, tolerance); + } + + /** + * 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 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 geom, double tolerance) { + return computeDistancePoints(geom, tolerance, -1.0); + } + + 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; + + if (geom.getDimension() == Dimension.P) { + return computeForPoints(geom, maxDistanceLimit); + } + //TODO: handle mixed geoms with points + Coordinate[] maxDistPtsEdge = computeForEdges(geom, tolerance, maxDistanceLimit); + + if (isBeyondLimit(distance(maxDistPtsEdge), maxDistanceLimit)) { + return maxDistPtsEdge; + } + + /** + * Polygonal query geometry may have an interior point as the farthest point. + */ + if (geom.getDimension() == Dimension.A) { + Coordinate[] maxDistPtsInterior = computeForAreaInterior(geom, tolerance); + if (maxDistPtsInterior != null + && distance(maxDistPtsInterior) > distance(maxDistPtsEdge)) { + return maxDistPtsInterior; + } + } + return maxDistPtsEdge; + } + + private Coordinate[] computeForPoints(Geometry geom, double maxDistanceLimit) { + double maxDist = -1.0;; + Coordinate[] maxDistPtsAB = null; + Iterator geomi = new GeometryCollectionIterator(geom); + while (geomi.hasNext()) { + Geometry geomElem = (Geometry) geomi.next(); + if (! (geomElem instanceof Point)) + continue; + + Coordinate pA = geomElem.getCoordinate(); + Coordinate pB = targetDistance.nearestPoint(pA); + double dist = pA.distance(pB); + + boolean isInterior = dist > 0 && targetDistance.isInterior(pA); + //-- check for interior point + if (isInterior) { + dist = 0; + pB = pA; + } + if (dist > maxDist) { + maxDist = dist; + maxDistPtsAB = pair(pA, pB); + } + if (isValidLimit(maxDistanceLimit) + && isBeyondLimit(maxDist, maxDistanceLimit)) { + break; + } + } + return maxDistPtsAB; + } + + private Coordinate[] computeForEdges(Geometry geom, double tolerance, double maxDistanceLimit) { + PriorityQueue segQueue = createSegQueue(geom); + + DHDSegment segMaxDist = null; + long iter = 0; + while (! segQueue.isEmpty()) { + 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.getMaxDistanceBound(); + double maxDist = segMaxBound.getMaxDistance(); + System.out.format("%s len: %f bound: %f maxDist: %f\n", + 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 + * 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 + */ + if (isValidLimit(maxDistanceLimit)) { + if (isWithinLimit(segMaxBound.getMaxDistanceBound(), maxDistanceLimit) + || isBeyondLimit(segMaxBound.getMaxDistance(), maxDistanceLimit) + ) { + break; + } + } + + /** + * Check for equal or collinear segments. + * If so, don't bisect the segment further. + * This greatly improves performance when the inputs + * have identical or collinear segments + * (in particular, the case when the inputs are identical). + */ + if (segMaxBound.getMaxDistance() == 0.0) { + if (isSameOrCollinear(segMaxBound)) + continue; + } + + //System.out.println(segMaxDist); + + /** + * If segment is longer than tolerance + * it might provide a better max distance point, + * so bisect and keep searching. + */ + if (tolerance > 0 + && segMaxBound.getLength() > tolerance) { + DHDSegment[] bisects = segMaxBound.bisect(targetDistance); + 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(); + + /** + * 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 = geom.getCoordinate(); + return pair(maxPt, maxPt); + } + + private boolean isSameOrCollinear(DHDSegment seg) { + CoordinateSequenceLocation f0 = targetDistance.nearestLocation(seg.p0); + CoordinateSequenceLocation f1 = targetDistance.nearestLocation(seg.p1); + 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; + } + + private static boolean isWithinLimit(double maxDist, double maxDistanceLimit) { + return maxDistanceLimit >= 0 && maxDist <= maxDistanceLimit; + } + + private void addNonInterior(DHDSegment segment, PriorityQueue segQueue) { + //-- discard segment if it is interior to a polygon + if (isInterior(segment)) { + return; + } + segQueue.add(segment); + } + + /** + * Tests if segment is fully in the interior of the target geometry polygons (if any). + * + * @param segment + * @return + */ + private boolean isInterior(DHDSegment segment) { + if (segment.getMaxDistance() > 0.0) { + return false; + } + return targetDistance.isInterior(segment.getEndpoint(0), segment.getEndpoint(1)); + } + + /** + * 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. + * + * 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 geom + * @param tolerance + * @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 geom, double tolerance) { + if (tolerance <= 0.0) + return null; + + //TODO: extract polygonal geoms from A + Geometry polygonal = geom; + + /** + * Optimization - skip if A interior cannot intersect B, + * and thus farther point must lie on A boundary + */ + if (polygonal.getEnvelopeInternal().disjoint(target.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. + */ + //TODO: add short-circuiting based on maxDistanceLimit? + + 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. + * This will cause the computed segment distance to be returned, + * which is preferred since it occurs on a vertex or edge. + */ + if (targetDistance.isInterior(ptA)) { + return null; + } + Coordinate ptB = targetDistance.nearestFacetPoint(ptA); + return pair(ptA, ptB); + } + + /** + * 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 geom + * @return the segment priority queue + */ + private PriorityQueue createSegQueue(Geometry geom) { + PriorityQueue priq = new PriorityQueue(); + geom.apply(new GeometryComponentFilter() { + + @Override + public void filter(Geometry geom) { + if (geom instanceof LineString) { + addSegments(geom.getCoordinates(), priq); + } + } + }); + return priq; + } + + /** + * Add segments to queue + * + * @param pts + * @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; + if (i == 0) { + seg = DHDSegment.create(pts[i], pts[i + 1], targetDistance); + } + else { + //-- avoiding recomputing prev pt distance + seg = DHDSegment.create(prevSeg, pts[i + 1], targetDistance); + } + prevSeg = seg; + + //-- don't add segment if it can't be further away then current max + if (segMaxDist == null + || seg.getMaxDistanceBound() > segMaxDist.getMaxDistance()) { + /** + * Don't add interior segments, since their distance must be zero. + */ + addNonInterior(seg, priq); + } + + if (segMaxDist == null + || seg.getMaxDistance() > segMaxDist.getMaxDistance()) { + segMaxDist = seg;; + } + //System.out.println(seg.distance()); + } + } + + private static class TargetDistance { + private IndexedFacetDistance distanceToFacets; + private boolean isArea; + private IndexedPointInPolygonsLocator ptInArea; + + public TargetDistance(Geometry geom) { + distanceToFacets = new IndexedFacetDistance(geom); + isArea = geom.getDimension() >= Dimension.A; + if (isArea) { + ptInArea = new IndexedPointInPolygonsLocator(geom); + } + } + + public CoordinateSequenceLocation nearestLocation(Coordinate p) { + return distanceToFacets.nearestLocation(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) { + 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 = distanceToFacets.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; + } + + } + + private static class DHDSegment implements Comparable { + + public static DHDSegment create(Coordinate p0, Coordinate p1, TargetDistance dist) { + DHDSegment seg = new DHDSegment(p0, p1); + seg.init(dist); + return seg; + } + + public static DHDSegment create(DHDSegment prevSeg, Coordinate p1, TargetDistance dist) { + DHDSegment seg = new DHDSegment(prevSeg.p1, p1); + seg.init(prevSeg.nearPt1, dist); + return seg; + } + + private Coordinate p0; + private Coordinate nearPt0; + private Coordinate p1; + private Coordinate nearPt1; + private double maxDistanceBound = Double.NEGATIVE_INFINITY; + private double maxDistance; + + 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; + computeMaxDistances(); + } + + private void init(TargetDistance dist) { + nearPt0 = dist.nearestPoint(p0); + nearPt1 = dist.nearestPoint(p1); + computeMaxDistances(); + } + + private void init(Coordinate nearest0, TargetDistance dist) { + nearPt0 = nearest0; + nearPt1 = dist.nearestPoint(p1); + computeMaxDistances(); + } + + public Coordinate getEndpoint(int index) { + return index == 0 ? p0 : p1; + } + + public double getLength() { + return p0.distance(p1); + } + + public double getMaxDistance() { + return maxDistance; + } + + public double getMaxDistanceBound() { + return maxDistanceBound; + } + + public Coordinate[] getMaxDistPts() { + double dist0 = p0.distance(nearPt0); + double dist1 = p1.distance(nearPt1); + if (dist0 > dist1) { + return pair(p0, nearPt0); + } + return pair(p1, nearPt1); + } + + /** + * Computes a least upper bound for the maximum distance to a segment. + */ + 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); + maxDistance = Math.max(dist0, dist1); + maxDistanceBound = maxDistance + getLength() / 2; + } + + public DHDSegment[] bisect(TargetDistance dist) { + Coordinate mid = new Coordinate( + (p0.x + p1.x) / 2, + (p0.y + p1.y) / 2 + ); + Coordinate nearPtMid = dist.nearestPoint(mid); + 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); + } + + public String toString() { + return WKTWriter.toLineString(p0, p1); + } + } +} 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/algorithm/construct/IndexedPointInPolygonsLocator.java b/modules/core/src/main/java/org/locationtech/jts/algorithm/locate/IndexedPointInPolygonsLocator.java similarity index 86% 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 ce99843041..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,12 +22,15 @@ /** * 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 * */ -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/operation/distance/CoordinateSequenceLocation.java b/modules/core/src/main/java/org/locationtech/jts/operation/distance/CoordinateSequenceLocation.java new file mode 100644 index 0000000000..d0d03c5a1b --- /dev/null +++ b/modules/core/src/main/java/org/locationtech/jts/operation/distance/CoordinateSequenceLocation.java @@ -0,0 +1,93 @@ +/* + * 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; + +/** + * A location on a {@link FacetSequence}. + * + * Location indexes are always the index of a sequence segment. + * 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 index of the final endpoint is + * normalized to 0. + * + * @author mdavis + * + */ +public class CoordinateSequenceLocation { + private CoordinateSequence seq; + private int index; + private Coordinate pt; + + public CoordinateSequenceLocation(CoordinateSequence seq, int index, Coordinate pt) { + this.seq = seq; + this.pt = pt; + this.index = index; + if (index >= seq.size()) + this.index = seq.size() - 1; + } + + public Coordinate getCoordinate() { + return pt; + } + + public boolean isSameSegment(CoordinateSequenceLocation f) { + 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() { + 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..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 @@ -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 CoordinateSequenceLocation nearestLocation(Coordinate p) + { + if (isPoint()) { + return new CoordinateSequenceLocation(pts, 0, pts.getCoordinate(0)); + } + return nearestLocationOnLine(p); + } + + private CoordinateSequenceLocation 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 2nd endpoint belongs to next segment + //-- except for last segment on non-closed sequence + if (dist == 0.0 && pt.equals2D(q1)) { + if (index < pts.size() - 1) { + index++; + } + //-- normalize index for a ring + index = normalize(pts, index); + } + if (minDistance <= 0.0) + break; + } + } + return new CoordinateSequenceLocation(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 e212ef4785..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 @@ -13,17 +13,19 @@ 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; /** - * 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 @@ -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); @@ -115,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, @@ -139,11 +125,11 @@ private Object[] nearestFacets(Geometry g) { } /** - * Computes the nearest points on the base geometry - * and the given 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 + * @return the nearest points on the target and the argument geometry */ public Coordinate[] nearestPoints(Geometry g) { @@ -153,8 +139,63 @@ public Coordinate[] nearestPoints(Geometry g) return fs1.nearestLocations(fs2); } + 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); + FacetSequence fsN = (FacetSequence) nearest; + return fsN.nearestLocation(p); + + } + /** - * Tests whether the base geometry lies within + * Computes the nearest point on the target geometry + * to a point. + * + * @param p the point coordinate + * @return the nearest point on the target geometry + */ + public Coordinate nearestPoint(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).getCoordinate(); + } + + /** + * Computes the distance from the target 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)); + } + + 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; + Coordinate[] loc = fs1.nearestLocations(fs2); + return loc[0].distance(loc[1]); + } + + /** + * Tests whether the target geometry lies within * a specified distance of the given geometry. * * @param g the geometry to test @@ -181,6 +222,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..2e83ef0ff3 --- /dev/null +++ b/modules/core/src/test/java/org/locationtech/jts/algorithm/distance/DirectedHausdorffDistanceTest.java @@ -0,0 +1,453 @@ +/* + * 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); } + + //-- empty inputs + + 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))"); + } + + //-------------------------------------- + //-- 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)", + "LINESTRING (0 0, 1 1)"); + } + + 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, "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); + } + + public void testPointsPolygon() + { + checkDistance("MULTIPOINT ((4 3), (2 8), (8 5))", "POLYGON ((6 9, 6 4, 9 1, 1 1, 6 9))", + "LINESTRING (2 8, 4.426966292134832 6.48314606741573)"); + } + + public void testLineSegments() + { + checkHausdorff("LINESTRING (0 0, 2 0)", "LINESTRING (0 0, 2 1)", + "LINESTRING (2 0, 2 1)"); + } + + public void testLineSegments2() + { + checkHausdorff("LINESTRING (0 0, 2 0)", "LINESTRING (0 1, 1 2, 2 1)", + "LINESTRING (1 0, 1 2)"); + } + + public void testLinePoints() + { + checkHausdorff("LINESTRING (0 0, 2 0)", "MULTIPOINT (0 2, 1 0, 2 1)", + "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.0); + } + + public void testLinesPolygon() + { + 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)"); + } + + 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, "LINESTRING (9 8, 6.3 7.1)"); + checkHausdorff(a, b, "LINESTRING (2 3, 5.5 3)"); + } + + public void testPolygonLineCrossingBoundaryResult() + { + checkDistance("POLYGON ((2 8, 8 2, 2 1, 2 8))", + "LINESTRING (6 5, 4 7, 0 0, 8 4)", + "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)", + "LINESTRING (4.555 2.989, 4.828 0.536)", 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, "LINESTRING (20 0, 17 3)"); + checkDistance(a, b, "LINESTRING (6.6796875 18, 11 19)"); + checkHausdorff(a, b, "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, "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, "LINESTRING (1 1, 5.4 3.2)"); + checkDistanceStartPtLen(b, a, + "LINESTRING (2.669921875 12.556640625, 5.446115154109589 13.818546660958905)", + 0.01); + } + + // 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, "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)"; + + checkHausdorff(wkt1, wkt2, "LINESTRING (6.53857421875 6.5382080078125, 6.53857421875 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, "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, "LINESTRING (3 5, 1 5)"); + checkDistance(wkt2, wkt1, "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); + } + + /** + * 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, 0.0); + } + + //----------------------------------------------------- + + 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() + { + String a = "MULTIPOINT ((1 9), (9 1))"; + String b = "MULTIPOINT ((1 1), (9 9))"; + 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, 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, 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, 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, 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, false); + checkFullyWithinDistance(a, b, 2, true); + checkFullyWithinDistance(a, b, 3, true); + } + + //====================================================================== + + private static final double TOLERANCE = 0.001; + + private void checkHausdorff(String wkt1, String wkt2, String wktExpected) + { + Geometry g1 = read(wkt1); + Geometry g2 = read(wkt2); + + Coordinate[] pts = DirectedHausdorffDistance.hausdorffDistancePoints(g1, g2); + Geometry result = g1.getFactory().createLineString(pts); + Geometry expected = read(wktExpected); + checkEqualExact(expected, result, TOLERANCE); + } + + private void checkDistance(String wkt1, String wkt2, String wktExpected) { + Geometry g1 = read(wkt1); + Geometry g2 = read(wkt2); + + Coordinate[] pts = DirectedHausdorffDistance.distancePoints(g1, g2); + Geometry result = g1.getFactory().createLineString(pts); + Geometry expected = read(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); + Geometry g2 = read(wkt2); + + Coordinate[] pts = DirectedHausdorffDistance.distancePoints(g1, g2); + Geometry result = g1.getFactory().createLineString(pts); + 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); + + double distResult = DirectedHausdorffDistance.distance(g1, g2, tolerance); + assertEquals(expectedDistance, distResult, TOLERANCE); + } + + private void checkDistance(String wkt1, String wkt2, + double expectedDistance) { + Geometry g1 = read(wkt1); + Geometry g2 = read(wkt2); + + double distResult = DirectedHausdorffDistance.distance(g1, g2); + assertEquals(expectedDistance, distResult, TOLERANCE); + } + + 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); + 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)); + } +} 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..d1d0876d5d --- /dev/null +++ b/modules/core/src/test/java/test/jts/perf/operation/distance/HausdorffDistancePerfTest.java @@ -0,0 +1,137 @@ +/* + * 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.DirectedHausdorffDistance; +import org.locationtech.jts.algorithm.distance.DiscreteHausdorffDistance; +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; + +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); + + //if (iter > 10) break; + } + } + + private void checkHausdorff(Geometry g, Geometry b) { + double distDHD = DiscreteHausdorffDistance.distance(g, b, 0.01); + + double distHD = DirectedHausdorffDistance.hausdorffDistance(g, b); + //-- 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.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); + 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; + } +} 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..3b7fb1693b --- /dev/null +++ b/modules/core/src/test/java/test/jts/perf/operation/distance/HausdorffDistanceStressTest.java @@ -0,0 +1,112 @@ +/* + * 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); + checkFullyWithinDistance(g, b); + + //if (iter > 10) break; + } + } + + private void checkHausdorff(Geometry g, Geometry b) { + double distDHD = DiscreteHausdorffDistance.distance(g, b, 0.01); + + double distHD = DirectedHausdorffDistance.hausdorffDistance(g, b); + //-- 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 void checkFullyWithinDistance(Geometry g, Geometry b) { + 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); + } + } + + 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; + } +}