[geos-commits] [SCM] GEOS branch main updated. 0d790d12b25a8e1d7310a55ec701e6bf54f73fb2
git at osgeo.org
git at osgeo.org
Mon Oct 5 09:35:41 PDT 2026
This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "GEOS".
The branch, main has been updated
via 0d790d12b25a8e1d7310a55ec701e6bf54f73fb2 (commit)
from c6af5613f1b5da92205e3cd45057499b90bfa0b8 (commit)
Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.
- Log -----------------------------------------------------------------
commit 0d790d12b25a8e1d7310a55ec701e6bf54f73fb2
Author: Jeroen Bloemscheer <jbloemscheer at gmail.com>
Date: Mon Oct 5 18:35:13 2026 +0200
Port JTS DirectedHausdorffDistance (#1514)
Co-authored-by: Dan Baston <dbaston at gmail.com>
diff --git a/NEWS.md b/NEWS.md
index 6891059d3..a6510e8d4 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,6 +1,10 @@
## Changes in 3.16.0
2027-xx-xx
+- New things:
+ - Add DirectedHausdorffDistance and C API GEOSDirectedHausdorffDistance /
+ GEOSSymmetricHausdorffDistance (Martin Davis, Jeroen Bloemscheer, Dan Baston)
+
- Fixes/Improvements:
- Curve overlay robustness improvements (GH-1513, Dan Baston)
- GeometrySplitter: Allow splitting a MultiSurface (GH-1520, Dan Baston)
diff --git a/Version.txt b/Version.txt
index c1bb857ba..50e4eb502 100644
--- a/Version.txt
+++ b/Version.txt
@@ -15,9 +15,9 @@ GEOS_PATCH_WORD=dev
# - Deleting interfaces / compatibility issues - bump CURRENT, others to zero
# ( THIS MUST BE CAREFULLY AVOIDED )
#
-CAPI_INTERFACE_CURRENT=23
+CAPI_INTERFACE_CURRENT=24
CAPI_INTERFACE_REVISION=0
-CAPI_INTERFACE_AGE=22
+CAPI_INTERFACE_AGE=23
# JTS Port
JTS_PORT=1.20.0
diff --git a/capi/geos_c.cpp b/capi/geos_c.cpp
index b2918e466..a5459ee92 100644
--- a/capi/geos_c.cpp
+++ b/capi/geos_c.cpp
@@ -409,6 +409,30 @@ extern "C" {
return GEOSHausdorffDistanceDensifyWithPoints_r(handle, g1, g2, densifyFrac, dist, p1x, p1y, p2x, p2y);
}
+ int
+ GEOSDirectedHausdorffDistance(const Geometry* g1, const Geometry* g2, double* dist)
+ {
+ return GEOSDirectedHausdorffDistance_r(handle, g1, g2, dist);
+ }
+
+ int
+ GEOSDirectedHausdorffDistanceWithPoints(const Geometry* g1, const Geometry* g2, double* dist, double* p1x, double* p1y, double* p2x, double* p2y)
+ {
+ return GEOSDirectedHausdorffDistanceWithPoints_r(handle, g1, g2, dist, p1x, p1y, p2x, p2y);
+ }
+
+ int
+ GEOSDirectedHausdorffDistanceWithin(const Geometry* g1, const Geometry* g2, double maxDistance)
+ {
+ return GEOSDirectedHausdorffDistanceWithin_r(handle, g1, g2, maxDistance);
+ }
+
+ int
+ GEOSSymmetricHausdorffDistance(const Geometry* g1, const Geometry* g2, double* dist)
+ {
+ return GEOSSymmetricHausdorffDistance_r(handle, g1, g2, dist);
+ }
+
int
GEOSFrechetDistance(const Geometry* g1, const Geometry* g2, double* dist)
{
diff --git a/capi/geos_c.h.in b/capi/geos_c.h.in
index 8121af536..00ec6b3ee 100644
--- a/capi/geos_c.h.in
+++ b/capi/geos_c.h.in
@@ -2184,6 +2184,36 @@ extern int GEOS_DLL GEOSHausdorffDistanceDensifyWithPoints_r(
double *p1x, double *p1y,
double *p2x, double *p2y);
+/** \see GEOSDirectedHausdorffDistance */
+extern int GEOS_DLL GEOSDirectedHausdorffDistance_r(
+ GEOSContextHandle_t handle,
+ const GEOSGeometry *g1,
+ const GEOSGeometry *g2,
+ double *dist);
+
+/** \see GEOSDirectedHausdorffDistanceWithPoints */
+extern int GEOS_DLL GEOSDirectedHausdorffDistanceWithPoints_r(
+ GEOSContextHandle_t handle,
+ const GEOSGeometry *g1,
+ const GEOSGeometry *g2,
+ double *dist,
+ double *p1x, double *p1y,
+ double *p2x, double *p2y);
+
+/** \see GEOSDirectedHausdorffDistanceWithin */
+extern int GEOS_DLL GEOSDirectedHausdorffDistanceWithin_r(
+ GEOSContextHandle_t handle,
+ const GEOSGeometry *g1,
+ const GEOSGeometry *g2,
+ double maxDistance);
+
+/** \see GEOSSymmetricHausdorffDistance */
+extern int GEOS_DLL GEOSSymmetricHausdorffDistance_r(
+ GEOSContextHandle_t handle,
+ const GEOSGeometry *g1,
+ const GEOSGeometry *g2,
+ double *dist);
+
/** \see GEOSFrechetDistance */
extern int GEOS_DLL GEOSFrechetDistance_r(
GEOSContextHandle_t handle,
@@ -4223,6 +4253,92 @@ extern int GEOS_DLL GEOSHausdorffDistanceDensifyWithPoints(
double* p1x, double* p1y,
double* p2x, double* p2y);
+/**
+* Calculate the directed Hausdorff distance from g1 to g2,
+* the largest distance from any point of g1 to g2.
+* Empty inputs write NaN.
+*
+* Approximated using a tolerance of the query envelope diameter / 1e4.
+*
+* @INPUT_CURVES_CONVERTED_TO_LINES@
+*
+* \param[in] g1 Query geometry
+* \param[in] g2 Target geometry
+* \param[out] dist Pointer to be filled in with distance result
+* \return 1 on success, 0 on exception.
+* \since 3.16
+*/
+extern int GEOS_DLL GEOSDirectedHausdorffDistance(
+ const GEOSGeometry *g1,
+ const GEOSGeometry *g2,
+ double *dist);
+
+/**
+* Calculate the directed Hausdorff distance from g1 to g2,
+* the largest distance from any point of g1 to g2,
+* and the coordinates of a point pair that realizes the distance
+* (a point of g1 and the nearest point of g2).
+* Empty inputs write NaN.
+*
+* Approximated using a tolerance of the query envelope diameter / 1e4.
+*
+* @INPUT_CURVES_CONVERTED_TO_LINES@
+*
+* \param[in] g1 Query geometry
+* \param[in] g2 Target geometry
+* \param[out] dist Pointer to be filled in with distance result
+* \param[out] p1x X coordinate of point on g1
+* \param[out] p1y Y coordinate of point on g1
+* \param[out] p2x X coordinate of point on g2
+* \param[out] p2y Y coordinate of point on g2
+* \return 1 on success, 0 on exception.
+* \since 3.16
+*/
+extern int GEOS_DLL GEOSDirectedHausdorffDistanceWithPoints(
+ const GEOSGeometry *g1,
+ const GEOSGeometry *g2,
+ double *dist,
+ double *p1x, double *p1y,
+ double *p2x, double *p2y);
+
+/**
+* Test whether every point of g1 lies within maxDistance of g2.
+* Empty inputs return 0.
+* Uses the same automatic approximation tolerance as GEOSDirectedHausdorffDistance.
+*
+* @INPUT_CURVES_CONVERTED_TO_LINES@
+*
+* \param[in] g1 Query geometry
+* \param[in] g2 Target geometry
+* \param[in] maxDistance the distance limit
+* \return 1 if fully within, 0 if not, 2 on exception.
+* \since 3.16
+*/
+extern int GEOS_DLL GEOSDirectedHausdorffDistanceWithin(
+ const GEOSGeometry *g1,
+ const GEOSGeometry *g2,
+ double maxDistance);
+
+/**
+* Calculate the Hausdorff distance between g1 and g2,
+* the larger of the two directed Hausdorff distances.
+* Empty inputs write NaN.
+*
+* Approximated using a tolerance of the query envelope diameter / 1e4.
+*
+* @INPUT_CURVES_CONVERTED_TO_LINES@
+*
+* \param[in] g1 Input geometry
+* \param[in] g2 Input geometry
+* \param[out] dist Pointer to be filled in with distance result
+* \return 1 on success, 0 on exception.
+* \since 3.16
+*/
+extern int GEOS_DLL GEOSSymmetricHausdorffDistance(
+ const GEOSGeometry *g1,
+ const GEOSGeometry *g2,
+ double *dist);
+
/**
* Calculate the
* [Frechet distance](https://en.wikipedia.org/wiki/Fr%C3%A9chet_distance)
diff --git a/capi/geos_ts_c.cpp b/capi/geos_ts_c.cpp
index a69e944ea..56218cb3b 100644
--- a/capi/geos_ts_c.cpp
+++ b/capi/geos_ts_c.cpp
@@ -27,6 +27,7 @@
#include <geos/algorithm/construct/MaximumInscribedCircle.h>
#include <geos/algorithm/construct/LargestEmptyCircle.h>
#include <geos/algorithm/distance/DiscreteHausdorffDistance.h>
+#include <geos/algorithm/distance/DirectedHausdorffDistance.h>
#include <geos/algorithm/distance/DiscreteFrechetDistance.h>
#include <geos/algorithm/hull/ConcaveHull.h>
#include <geos/algorithm/hull/ConcaveHullOfPolygons.h>
@@ -135,6 +136,7 @@
#include <cstring>
#include <fstream>
#include <iostream>
+#include <limits>
#include <optional>
#include <sstream>
#include <string>
@@ -232,6 +234,7 @@ using geos::io::GeoJSONWriter;
using geos::algorithm::distance::DiscreteFrechetDistance;
using geos::algorithm::distance::DiscreteHausdorffDistance;
+using geos::algorithm::distance::DirectedHausdorffDistance;
using geos::algorithm::hull::ConcaveHull;
using geos::algorithm::hull::ConcaveHullOfPolygons;
@@ -1287,6 +1290,66 @@ extern "C" {
});
}
+ int
+ GEOSDirectedHausdorffDistance_r(GEOSContextHandle_t extHandle, const Geometry* g1, const Geometry* g2, double* dist)
+ {
+ return execute(extHandle, 0, [&]() {
+ const auto input1 = convertToLineIfNeeded(extHandle, g1);
+ const auto input2 = convertToLineIfNeeded(extHandle, g2);
+
+ *dist = DirectedHausdorffDistance::distance(*input1, *input2);
+ return 1;
+ });
+ }
+
+ int
+ GEOSDirectedHausdorffDistanceWithPoints_r(GEOSContextHandle_t extHandle,
+ const Geometry* g1, const Geometry* g2,
+ double* dist, double* p1x, double* p1y,
+ double* p2x, double* p2y)
+ {
+ return execute(extHandle, 0, [&]() {
+ const auto input1 = convertToLineIfNeeded(extHandle, g1);
+ const auto input2 = convertToLineIfNeeded(extHandle, g2);
+
+ auto pts = DirectedHausdorffDistance::distancePoints(*input1, *input2);
+ if (!pts) {
+ *dist = std::numeric_limits<double>::quiet_NaN();
+ *p1x = *p1y = *p2x = *p2y = std::numeric_limits<double>::quiet_NaN();
+ return 1;
+ }
+ *dist = (*pts)[0].distance((*pts)[1]);
+ *p1x = (*pts)[0].x;
+ *p1y = (*pts)[0].y;
+ *p2x = (*pts)[1].x;
+ *p2y = (*pts)[1].y;
+ return 1;
+ });
+ }
+
+ int
+ GEOSDirectedHausdorffDistanceWithin_r(GEOSContextHandle_t extHandle,
+ const Geometry* g1, const Geometry* g2, double maxDistance)
+ {
+ return execute(extHandle, 2, [&]() {
+ const auto input1 = convertToLineIfNeeded(extHandle, g1);
+ const auto input2 = convertToLineIfNeeded(extHandle, g2);
+ return DirectedHausdorffDistance::isFullyWithinDistance(*input1, *input2, maxDistance) ? 1 : 0;
+ });
+ }
+
+ int
+ GEOSSymmetricHausdorffDistance_r(GEOSContextHandle_t extHandle, const Geometry* g1, const Geometry* g2, double* dist)
+ {
+ return execute(extHandle, 0, [&]() {
+ const auto input1 = convertToLineIfNeeded(extHandle, g1);
+ const auto input2 = convertToLineIfNeeded(extHandle, g2);
+
+ *dist = DirectedHausdorffDistance::hausdorffDistance(*input1, *input2);
+ return 1;
+ });
+ }
+
int
GEOSFrechetDistance_r(GEOSContextHandle_t extHandle, const Geometry* g1, const Geometry* g2, double* dist)
{
diff --git a/include/geos/algorithm/distance/DirectedHausdorffDistance.h b/include/geos/algorithm/distance/DirectedHausdorffDistance.h
new file mode 100644
index 000000000..8ccbb4a22
--- /dev/null
+++ b/include/geos/algorithm/distance/DirectedHausdorffDistance.h
@@ -0,0 +1,300 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2026 Martin Davis
+ * Copyright (C) 2026 Jeroen Bloemscheer
+ *
+ * This is free software; you can redistribute and/or modify it under
+ * the terms of the GNU Lesser General Public Licence as published
+ * by the Free Software Foundation.
+ * See the COPYING file for more information.
+ *
+ **********************************************************************
+ *
+ * Last port: algorithm/distance/DirectedHausdorffDistance.java (aff11591)
+ *
+ **********************************************************************/
+
+#pragma once
+
+#include <geos/export.h>
+#include <geos/geom/Coordinate.h>
+
+#include <array>
+#include <memory>
+#include <optional>
+
+#ifdef _MSC_VER
+#pragma warning(push)
+#pragma warning(disable: 4251) // warning C4251: needs to have dll-interface to be used by clients of class
+#endif
+
+namespace geos {
+namespace geom {
+class Envelope;
+class Geometry;
+}
+}
+
+namespace geos {
+namespace algorithm {
+namespace distance {
+
+/**
+ * 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) = max a ∈ A (max b ∈ 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 a symmetric distance metric defined as:
+ *
+ * HD(A,B) = max(DHD(A,B), DHD(B,A))
+ *
+ * This can be computed via the
+ * 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 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 farthestPoints(Geometry)
+ * or 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 DiscreteHausdorffDistance.
+ *
+ * \author Martin Davis
+ */
+class GEOS_DLL DirectedHausdorffDistance {
+public:
+ /// JTS returns Coordinate[]; std::array is the C++ equivalent pair.
+ using PointPair = std::array<geom::CoordinateXY, 2>;
+
+ /**
+ * 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
+ */
+ static double distance(const geom::Geometry& a, const geom::Geometry& 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 accuracy distance tolerance
+ * @return the directed Hausdorff distance,
+ * or NaN if an input is empty
+ */
+ static double distance(const geom::Geometry& a, const geom::Geometry& b,
+ double 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
+ * @return a pair of points [ptA, ptB] demonstrating the distance,
+ * or empty if an input is empty
+ */
+ static std::optional<PointPair> distancePoints(
+ const geom::Geometry& a, const geom::Geometry& b);
+
+ /**
+ * 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 empty if an input is empty
+ */
+ static std::optional<PointPair> distancePoints(
+ const geom::Geometry& a, const geom::Geometry& b, double tolerance);
+
+ /**
+ * Computes the symmetric Hausdorff distance between two geometries.
+ * This is 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
+ */
+ static double hausdorffDistance(const geom::Geometry& a, const geom::Geometry& b);
+
+ /**
+ * Computes a pair of points which attain the symmetric Hausdorff distance
+ * between two geometries.
+ * This is 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 empty if an input is empty
+ */
+ static std::optional<PointPair> hausdorffDistancePoints(
+ const geom::Geometry& a, const geom::Geometry& b);
+
+ /**
+ * Computes whether a query geometry lies fully within a given 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
+ */
+ static bool isFullyWithinDistance(
+ const geom::Geometry& a, const geom::Geometry& b, double maxDistance);
+
+ /**
+ * Computes whether a query geometry lies fully within a given 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
+ */
+ static bool isFullyWithinDistance(
+ const geom::Geometry& a, const geom::Geometry& b,
+ double maxDistance, double tolerance);
+
+ /**
+ * Create a new instance for a target geometry.
+ *
+ * @param geom the geometry to compute the distance from
+ */
+ explicit DirectedHausdorffDistance(const geom::Geometry& geom);
+
+ DirectedHausdorffDistance(const DirectedHausdorffDistance&) = delete;
+ DirectedHausdorffDistance& operator=(const DirectedHausdorffDistance&) = delete;
+
+ ~DirectedHausdorffDistance();
+
+ /**
+ * 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 empty.
+ *
+ * @param geom the query geometry
+ * @return a pair of points [ptA, ptB] attaining the distance,
+ * or empty if an input is empty
+ */
+ std::optional<PointPair> farthestPoints(const geom::Geometry& geom) const;
+
+ /**
+ * 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 empty.
+ *
+ * @param geom the query geometry
+ * @param tolerance the approximation distance tolerance
+ * @return a pair of points [ptA, ptB] attaining the distance,
+ * or empty if an input is empty
+ */
+ std::optional<PointPair> farthestPoints(const geom::Geometry& geom, double tolerance) const;
+
+ /**
+ * Tests whether a query geometry lies fully within a given 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
+ */
+ bool isFullyWithinDistance(const geom::Geometry& geom, double maxDistance) const;
+
+ /**
+ * Tests whether a query geometry lies fully within a given 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
+ */
+ bool isFullyWithinDistance(
+ const geom::Geometry& geom, double maxDistance, double tolerance) const;
+
+private:
+ class TargetDistance;
+ friend class DHDSegment;
+
+ static double pairDistance(const std::optional<PointPair>& pts);
+ static PointPair pair(const geom::CoordinateXY& p0, const geom::CoordinateXY& p1);
+ static double computeTolerance(const geom::Geometry& geom);
+ static bool isBeyond(
+ const geom::Envelope& envA, const geom::Envelope& envB, double maxDistance);
+ static bool isValidLimit(double limit);
+ static bool isBeyondLimit(double maxDist, double maxDistanceLimit);
+ static bool isWithinLimit(double maxDist, double maxDistanceLimit);
+
+ std::optional<PointPair> computeDistancePoints(
+ const geom::Geometry& geom, double tolerance, double maxDistanceLimit) const;
+ std::optional<PointPair> computeForPoints(
+ const geom::Geometry& geom, double maxDistanceLimit) const;
+ std::optional<PointPair> computeForEdges(
+ const geom::Geometry& geom, double tolerance, double maxDistanceLimit) const;
+ std::optional<PointPair> computeForAreaInterior(
+ const geom::Geometry& geom, double tolerance) const;
+
+ const geom::Geometry& target;
+ std::unique_ptr<TargetDistance> targetDistance;
+};
+
+} // namespace distance
+} // namespace algorithm
+} // namespace geos
+
+#ifdef _MSC_VER
+#pragma warning(pop)
+#endif
diff --git a/include/geos/operation/distance/CoordinateSequenceLocation.h b/include/geos/operation/distance/CoordinateSequenceLocation.h
new file mode 100644
index 000000000..8686961ac
--- /dev/null
+++ b/include/geos/operation/distance/CoordinateSequenceLocation.h
@@ -0,0 +1,76 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2026 Martin Davis
+ *
+ * This is free software; you can redistribute and/or modify it under
+ * the terms of the GNU Lesser General Public Licence as published
+ * by the Free Software Foundation.
+ * See the COPYING file for more information.
+ *
+ **********************************************************************
+ *
+ * Last port: operation/distance/CoordinateSequenceLocation.java (aff11591)
+ *
+ **********************************************************************/
+
+#pragma once
+
+#include <geos/export.h>
+#include <geos/geom/Coordinate.h>
+
+#include <cstddef>
+
+namespace geos {
+namespace geom {
+class CoordinateSequence;
+}
+}
+
+namespace geos {
+namespace operation {
+namespace distance {
+
+/**
+ * A location on a 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 Martin Davis
+ */
+class GEOS_DLL CoordinateSequenceLocation {
+public:
+
+ CoordinateSequenceLocation(const geom::CoordinateSequence* seq,
+ std::size_t index,
+ const geom::CoordinateXY& pt);
+
+ const geom::CoordinateXY& getCoordinate() const;
+
+ std::size_t getIndex() const;
+
+ bool isSameSegment(const CoordinateSequenceLocation& f) const;
+
+ const geom::CoordinateXY& getEndPoint(int i) const;
+
+ std::size_t normalize(std::size_t index) const;
+
+private:
+
+ const geom::CoordinateSequence* seq;
+ std::size_t index;
+ geom::CoordinateXY pt;
+
+ bool isNext(std::size_t index, std::size_t index1) const;
+};
+
+} // namespace distance
+} // namespace operation
+} // namespace geos
diff --git a/include/geos/operation/distance/FacetSequence.h b/include/geos/operation/distance/FacetSequence.h
index 61f9fd65b..2f09acd46 100644
--- a/include/geos/operation/distance/FacetSequence.h
+++ b/include/geos/operation/distance/FacetSequence.h
@@ -12,7 +12,7 @@
*
**********************************************************************
*
- * Last port: operation/distance/FacetSequence.java (f6187ee2 JTS-1.14)
+ * Last port: operation/distance/FacetSequence.java (aff11591)
*
**********************************************************************/
@@ -20,6 +20,9 @@
#include <geos/geom/Envelope.h>
#include <geos/geom/Coordinate.h>
+#include <geos/operation/distance/CoordinateSequenceLocation.h>
+
+#include <vector>
namespace geos {
namespace geom {
@@ -61,6 +64,10 @@ private:
const geom::Coordinate& q0, const geom::Coordinate &q1,
std::vector<geom::Coordinate> *locs) const;
+ CoordinateSequenceLocation nearestLocationOnLine(const geom::CoordinateXY& pt) const;
+
+ static std::size_t normalize(const geom::CoordinateSequence& pts, std::size_t index);
+
void computeEnvelope();
public:
@@ -82,13 +89,21 @@ public:
* and another sequence.
* The locations are presented in the same order as the input sequences.
*
- * @return a pair of {@link Coordinate}s for the nearest points
+ * @return a pair of Coordinates for the nearest points
*/
std::vector<geom::Coordinate> nearestLocations(const FacetSequence& facetSeq) const;
+ /**
+ * Computes the location of the nearest point of this sequence
+ * to a query point.
+ *
+ * @param p the query point
+ * @return the nearest location on this sequence
+ */
+ CoordinateSequenceLocation nearestLocation(const geom::CoordinateXY& p) const;
+
};
}
}
}
-
diff --git a/include/geos/operation/distance/IndexedFacetDistance.h b/include/geos/operation/distance/IndexedFacetDistance.h
index 6739efbf0..46e23c034 100644
--- a/include/geos/operation/distance/IndexedFacetDistance.h
+++ b/include/geos/operation/distance/IndexedFacetDistance.h
@@ -12,12 +12,13 @@
*
**********************************************************************
*
- * Last port: operation/distance/IndexedFacetDistance.java (f6187ee2 JTS-1.14)
+ * Last port: operation/distance/IndexedFacetDistance.java (aff11591)
*
**********************************************************************/
#pragma once
+#include <geos/operation/distance/CoordinateSequenceLocation.h>
#include <geos/operation/distance/FacetSequenceTreeBuilder.h>
namespace geos {
@@ -102,6 +103,21 @@ public:
/// \return the nearest points
std::unique_ptr<geom::CoordinateSequence> nearestPoints(const geom::Geometry* g) const;
+ /**
+ * Computes the nearest point on the target geometry
+ * to a point.
+ *
+ * @param p the point coordinate
+ * @return the nearest point on the target geometry
+ */
+ geom::CoordinateXY nearestPoint(const geom::CoordinateXY& p) const;
+
+ double distance(const geom::CoordinateXY& p) const;
+
+ double distance(const geom::CoordinateXY& p0, const geom::CoordinateXY& p1) const;
+
+ CoordinateSequenceLocation nearestLocation(const geom::CoordinateXY& p) const;
+
private:
struct FacetDistance {
diff --git a/src/algorithm/distance/DirectedHausdorffDistance.cpp b/src/algorithm/distance/DirectedHausdorffDistance.cpp
new file mode 100644
index 000000000..a5a4aa6ab
--- /dev/null
+++ b/src/algorithm/distance/DirectedHausdorffDistance.cpp
@@ -0,0 +1,659 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2026 Martin Davis
+ * Copyright (C) 2026 Jeroen Bloemscheer
+ *
+ * This is free software; you can redistribute and/or modify it under
+ * the terms of the GNU Lesser General Public Licence as published
+ * by the Free Software Foundation.
+ * See the COPYING file for more information.
+ *
+ **********************************************************************
+ *
+ * Last port: algorithm/distance/DirectedHausdorffDistance.java (aff11591)
+ *
+ **********************************************************************/
+
+#include <geos/algorithm/distance/DirectedHausdorffDistance.h>
+
+#include <geos/algorithm/construct/IndexedPointInPolygonsLocator.h>
+#include <geos/algorithm/construct/LargestEmptyCircle.h>
+#include <geos/geom/Coordinate.h>
+#include <geos/geom/CoordinateSequence.h>
+#include <geos/geom/Dimension.h>
+#include <geos/geom/Envelope.h>
+#include <geos/geom/Geometry.h>
+#include <geos/geom/GeometryComponentFilter.h>
+#include <geos/geom/LineString.h>
+#include <geos/geom/Location.h>
+#include <geos/geom/Point.h>
+#include <geos/operation/distance/IndexedFacetDistance.h>
+#include <geos/util/IllegalArgumentException.h>
+
+#include <cmath>
+#include <limits>
+#include <queue>
+#include <utility>
+#include <vector>
+
+using geos::algorithm::construct::IndexedPointInPolygonsLocator;
+using geos::algorithm::construct::LargestEmptyCircle;
+using geos::geom::CoordinateXY;
+using geos::geom::Dimension;
+using geos::geom::Envelope;
+using geos::geom::Geometry;
+using geos::geom::GeometryComponentFilter;
+using geos::geom::LineString;
+using geos::geom::Location;
+using geos::geom::Point;
+using geos::operation::distance::IndexedFacetDistance;
+
+namespace geos {
+namespace algorithm {
+namespace distance {
+
+namespace {
+
+constexpr double EMPTY_DISTANCE = std::numeric_limits<double>::quiet_NaN();
+constexpr double AUTO_TOLERANCE_FACTOR = 1.0e4;
+constexpr double AREA_INTERIOR_TOLERANCE_FACTOR = 20;
+constexpr double FULLY_WITHIN_TOLERANCE_FACTOR = 10 * AUTO_TOLERANCE_FACTOR;
+
+} // namespace
+
+class DirectedHausdorffDistance::TargetDistance {
+public:
+ explicit TargetDistance(const Geometry& geom)
+ : distanceToFacets(&geom)
+ , isArea(geom.getDimension() >= Dimension::A)
+ {
+ if (isArea) {
+ ptInArea = std::make_unique<IndexedPointInPolygonsLocator>(geom);
+ }
+ }
+
+ CoordinateXY nearestFacetPoint(const CoordinateXY& p) const
+ {
+ CoordinateXY np = distanceToFacets.nearestPoint(p);
+ return np;
+ }
+
+ CoordinateXY nearestPoint(const CoordinateXY& p) const
+ {
+ if (ptInArea) {
+ if (ptInArea->locate(&p) != Location::EXTERIOR) {
+ return p;
+ }
+ }
+ return distanceToFacets.nearestPoint(p);
+ }
+
+ bool isInterior(const CoordinateXY& p) const
+ {
+ if (!isArea) {
+ return false;
+ }
+ return ptInArea->locate(&p) == Location::INTERIOR;
+ }
+
+ bool isInterior(const CoordinateXY& p0, const CoordinateXY& p1) const
+ {
+ if (!isArea) {
+ return false;
+ }
+ double segDist = distanceToFacets.distance(p0, p1);
+ if (segDist == 0.0) {
+ return false;
+ }
+ return isInterior(p0);
+ }
+
+ /// JTS isSameOrCollinear: both endpoints project onto the same target segment.
+ bool isSameOrCollinear(const CoordinateXY& p0, const CoordinateXY& p1) const
+ {
+ const auto f0 = distanceToFacets.nearestLocation(p0);
+ const auto f1 = distanceToFacets.nearestLocation(p1);
+ return f0.isSameSegment(f1);
+ }
+
+private:
+ IndexedFacetDistance distanceToFacets;
+ bool isArea;
+ std::unique_ptr<IndexedPointInPolygonsLocator> ptInArea;
+};
+
+class DHDSegment {
+public:
+ static DHDSegment create(const CoordinateXY& p0, const CoordinateXY& p1,
+ DirectedHausdorffDistance::TargetDistance& dist)
+ {
+ DHDSegment seg(p0, p1);
+ seg.init(dist);
+ return seg;
+ }
+
+ static DHDSegment create(const DHDSegment& prevSeg, const CoordinateXY& p1,
+ DirectedHausdorffDistance::TargetDistance& dist)
+ {
+ DHDSegment seg(prevSeg.p1, p1);
+ seg.init(prevSeg.nearPt1, dist);
+ return seg;
+ }
+
+ CoordinateXY getEndpoint(int index) const
+ {
+ return index == 0 ? p0 : p1;
+ }
+
+ double getLength() const
+ {
+ return p0.distance(p1);
+ }
+
+ double getMaxDistance() const
+ {
+ return maxDistance;
+ }
+
+ double getMaxDistanceBound() const
+ {
+ return maxDistanceBound;
+ }
+
+ DirectedHausdorffDistance::PointPair getMaxDistPts() const
+ {
+ double dist0 = p0.distance(nearPt0);
+ double dist1 = p1.distance(nearPt1);
+ if (dist0 > dist1) {
+ return DirectedHausdorffDistance::PointPair{p0, nearPt0};
+ }
+ return DirectedHausdorffDistance::PointPair{p1, nearPt1};
+ }
+
+ std::array<DHDSegment, 2> bisect(DirectedHausdorffDistance::TargetDistance& dist) const
+ {
+ CoordinateXY mid((p0.x + p1.x) / 2.0, (p0.y + p1.y) / 2.0);
+ CoordinateXY nearPtMid = dist.nearestPoint(mid);
+ return {
+ DHDSegment(p0, nearPt0, mid, nearPtMid),
+ DHDSegment(mid, nearPtMid, p1, nearPt1)
+ };
+ }
+
+ bool operator<(const DHDSegment& other) const
+ {
+ // priority_queue is a max-heap: larger bound first
+ return maxDistanceBound < other.maxDistanceBound;
+ }
+
+ CoordinateXY p0;
+ CoordinateXY nearPt0;
+ CoordinateXY p1;
+ CoordinateXY nearPt1;
+
+private:
+ DHDSegment(const CoordinateXY& p0_, const CoordinateXY& p1_)
+ : p0(p0_)
+ , p1(p1_)
+ , maxDistanceBound(-std::numeric_limits<double>::infinity())
+ , maxDistance(0.0)
+ {}
+
+ DHDSegment(const CoordinateXY& p0_, const CoordinateXY& nearPt0_,
+ const CoordinateXY& p1_, const CoordinateXY& nearPt1_)
+ : p0(p0_)
+ , nearPt0(nearPt0_)
+ , p1(p1_)
+ , nearPt1(nearPt1_)
+ , maxDistanceBound(-std::numeric_limits<double>::infinity())
+ , maxDistance(0.0)
+ {
+ computeMaxDistances();
+ }
+
+ void init(DirectedHausdorffDistance::TargetDistance& dist)
+ {
+ nearPt0 = dist.nearestPoint(p0);
+ nearPt1 = dist.nearestPoint(p1);
+ computeMaxDistances();
+ }
+
+ void init(const CoordinateXY& nearest0, DirectedHausdorffDistance::TargetDistance& dist)
+ {
+ nearPt0 = nearest0;
+ nearPt1 = dist.nearestPoint(p1);
+ computeMaxDistances();
+ }
+
+ void computeMaxDistances()
+ {
+ double dist0 = p0.distance(nearPt0);
+ double dist1 = p1.distance(nearPt1);
+ maxDistance = std::max(dist0, dist1);
+ maxDistanceBound = maxDistance + getLength() / 2.0;
+ }
+
+ double maxDistanceBound;
+ double maxDistance;
+};
+
+/* static */
+double
+DirectedHausdorffDistance::pairDistance(const std::optional<PointPair>& pts)
+{
+ if (!pts) {
+ return EMPTY_DISTANCE;
+ }
+ return (*pts)[0].distance((*pts)[1]);
+}
+
+/* static */
+DirectedHausdorffDistance::PointPair
+DirectedHausdorffDistance::pair(const CoordinateXY& p0, const CoordinateXY& p1)
+{
+ return PointPair{p0, p1};
+}
+
+/* static */
+double
+DirectedHausdorffDistance::computeTolerance(const Geometry& geom)
+{
+ return geom.getEnvelopeInternal()->getDiameter() / AUTO_TOLERANCE_FACTOR;
+}
+
+/* static */
+bool
+DirectedHausdorffDistance::isBeyond(
+ const Envelope& envA, const Envelope& envB, double maxDistance)
+{
+ if (envA.isNull() || envB.isNull()) {
+ return false;
+ }
+ return envA.getMinX() < envB.getMinX() - maxDistance
+ || envA.getMinY() < envB.getMinY() - maxDistance
+ || envA.getMaxX() > envB.getMaxX() + maxDistance
+ || envA.getMaxY() > envB.getMaxY() + maxDistance;
+}
+
+/* static */
+bool
+DirectedHausdorffDistance::isValidLimit(double limit)
+{
+ return limit >= 0.0;
+}
+
+/* static */
+bool
+DirectedHausdorffDistance::isBeyondLimit(double maxDist, double maxDistanceLimit)
+{
+ return maxDistanceLimit >= 0 && maxDist > maxDistanceLimit;
+}
+
+/* static */
+bool
+DirectedHausdorffDistance::isWithinLimit(double maxDist, double maxDistanceLimit)
+{
+ return maxDistanceLimit >= 0 && maxDist <= maxDistanceLimit;
+}
+
+/* static */
+double
+DirectedHausdorffDistance::distance(const Geometry& a, const Geometry& b)
+{
+ DirectedHausdorffDistance hd(b);
+ return pairDistance(hd.farthestPoints(a));
+}
+
+/* static */
+double
+DirectedHausdorffDistance::distance(const Geometry& a, const Geometry& b, double tolerance)
+{
+ DirectedHausdorffDistance hd(b);
+ return pairDistance(hd.farthestPoints(a, tolerance));
+}
+
+/* static */
+std::optional<DirectedHausdorffDistance::PointPair>
+DirectedHausdorffDistance::distancePoints(const Geometry& a, const Geometry& b)
+{
+ DirectedHausdorffDistance dhd(b);
+ return dhd.farthestPoints(a);
+}
+
+/* static */
+std::optional<DirectedHausdorffDistance::PointPair>
+DirectedHausdorffDistance::distancePoints(
+ const Geometry& a, const Geometry& b, double tolerance)
+{
+ DirectedHausdorffDistance dhd(b);
+ return dhd.farthestPoints(a, tolerance);
+}
+
+/* static */
+std::optional<DirectedHausdorffDistance::PointPair>
+DirectedHausdorffDistance::hausdorffDistancePoints(const Geometry& a, const Geometry& b)
+{
+ DirectedHausdorffDistance hdAB(b);
+ auto ptsAB = hdAB.farthestPoints(a);
+ DirectedHausdorffDistance hdBA(a);
+ auto ptsBA = hdBA.farthestPoints(b);
+
+ if (!ptsAB) {
+ return ptsBA ? std::optional<PointPair>(pair((*ptsBA)[1], (*ptsBA)[0])) : std::nullopt;
+ }
+ if (!ptsBA) {
+ return ptsAB;
+ }
+ if (pairDistance(ptsBA) > pairDistance(ptsAB)) {
+ return pair((*ptsBA)[1], (*ptsBA)[0]);
+ }
+ return ptsAB;
+}
+
+/* static */
+double
+DirectedHausdorffDistance::hausdorffDistance(const Geometry& a, const Geometry& b)
+{
+ return pairDistance(hausdorffDistancePoints(a, b));
+}
+
+/* static */
+bool
+DirectedHausdorffDistance::isFullyWithinDistance(
+ const Geometry& a, const Geometry& b, double maxDistance)
+{
+ DirectedHausdorffDistance hd(b);
+ return hd.isFullyWithinDistance(a, maxDistance);
+}
+
+/* static */
+bool
+DirectedHausdorffDistance::isFullyWithinDistance(
+ const Geometry& a, const Geometry& b, double maxDistance, double tolerance)
+{
+ DirectedHausdorffDistance hd(b);
+ return hd.isFullyWithinDistance(a, maxDistance, tolerance);
+}
+
+DirectedHausdorffDistance::DirectedHausdorffDistance(const Geometry& geom)
+ : target(geom)
+ , targetDistance(std::make_unique<TargetDistance>(geom))
+{}
+
+DirectedHausdorffDistance::~DirectedHausdorffDistance() = default;
+
+bool
+DirectedHausdorffDistance::isFullyWithinDistance(
+ const Geometry& geom, double maxDistance) const
+{
+ double tolerance = maxDistance / FULLY_WITHIN_TOLERANCE_FACTOR;
+ return isFullyWithinDistance(geom, maxDistance, tolerance);
+}
+
+bool
+DirectedHausdorffDistance::isFullyWithinDistance(
+ const Geometry& geom, double maxDistance, double tolerance) const
+{
+ if (geom.isEmpty() || target.isEmpty()) {
+ return false;
+ }
+ if (isBeyond(*geom.getEnvelopeInternal(), *target.getEnvelopeInternal(), maxDistance)) {
+ return false;
+ }
+ auto maxDistCoords = computeDistancePoints(geom, tolerance, maxDistance);
+ if (!maxDistCoords) {
+ return false;
+ }
+ return pairDistance(maxDistCoords) <= maxDistance;
+}
+
+std::optional<DirectedHausdorffDistance::PointPair>
+DirectedHausdorffDistance::farthestPoints(const Geometry& geom) const
+{
+ return farthestPoints(geom, computeTolerance(geom));
+}
+
+std::optional<DirectedHausdorffDistance::PointPair>
+DirectedHausdorffDistance::farthestPoints(const Geometry& geom, double tolerance) const
+{
+ return computeDistancePoints(geom, tolerance, -1.0);
+}
+
+std::optional<DirectedHausdorffDistance::PointPair>
+DirectedHausdorffDistance::computeDistancePoints(
+ const Geometry& geom, double tolerance, double maxDistanceLimit) const
+{
+ if (tolerance < 0.0) {
+ throw util::IllegalArgumentException("Tolerance must be non-negative");
+ }
+ if (geom.isEmpty() || target.isEmpty()) {
+ return std::nullopt;
+ }
+ if (geom.getDimension() == Dimension::P) {
+ return computeForPoints(geom, maxDistanceLimit);
+ }
+
+ auto maxDistPtsEdge = computeForEdges(geom, tolerance, maxDistanceLimit);
+ if (isBeyondLimit(pairDistance(maxDistPtsEdge), maxDistanceLimit)) {
+ return maxDistPtsEdge;
+ }
+ if (geom.getDimension() == Dimension::A) {
+ auto maxDistPtsInterior = computeForAreaInterior(geom, tolerance);
+ if (maxDistPtsInterior
+ && pairDistance(maxDistPtsInterior) > pairDistance(maxDistPtsEdge)) {
+ return maxDistPtsInterior;
+ }
+ }
+ return maxDistPtsEdge;
+}
+
+std::optional<DirectedHausdorffDistance::PointPair>
+DirectedHausdorffDistance::computeForPoints(
+ const Geometry& geom, double maxDistanceLimit) const
+{
+ double maxDist = -1.0;
+ std::optional<PointPair> maxDistPtsAB;
+
+ struct PointFilter : public GeometryComponentFilter {
+ PointFilter(TargetDistance& td, double& md, std::optional<PointPair>& pts,
+ double limit)
+ : targetDistance(td), maxDist(md), maxDistPtsAB(pts), maxDistanceLimit(limit), done(false)
+ {}
+
+ bool isDone() override {
+ return done;
+ }
+
+ void filter_ro(const Geometry* geomElem) override
+ {
+ if (geomElem->getGeometryTypeId() != geom::GEOS_POINT) {
+ return;
+ }
+
+ const CoordinateXY* pA = geomElem->getCoordinate();
+ CoordinateXY pB = targetDistance.nearestPoint(*pA);
+ double dist = pA->distance(pB);
+ bool interior = dist > 0 && targetDistance.isInterior(*pA);
+ if (interior) {
+ dist = 0;
+ pB = *pA;
+ }
+ if (dist > maxDist) {
+ maxDist = dist;
+ maxDistPtsAB = DirectedHausdorffDistance::pair(*pA, pB);
+ }
+ if (DirectedHausdorffDistance::isValidLimit(maxDistanceLimit)
+ && DirectedHausdorffDistance::isBeyondLimit(maxDist, maxDistanceLimit)) {
+ done = true;
+ }
+ }
+
+ TargetDistance& targetDistance;
+ double& maxDist;
+ std::optional<PointPair>& maxDistPtsAB;
+ double maxDistanceLimit;
+ bool done;
+ };
+
+ PointFilter filter(*targetDistance, maxDist, maxDistPtsAB, maxDistanceLimit);
+ geom.apply_ro(&filter);
+ return maxDistPtsAB;
+}
+
+std::optional<DirectedHausdorffDistance::PointPair>
+DirectedHausdorffDistance::computeForEdges(
+ const Geometry& geom, double tolerance, double maxDistanceLimit) const
+{
+ std::priority_queue<DHDSegment> segQueue;
+
+ struct LineFilter : public GeometryComponentFilter {
+ LineFilter(std::priority_queue<DHDSegment>& q, TargetDistance& td)
+ : segQueue(q), targetDistance(td)
+ {}
+
+ void filter_ro(const Geometry* g) override
+ {
+ const LineString* ls = dynamic_cast<const LineString*>(g);
+ if (!ls || ls->isEmpty()) {
+ return;
+ }
+ const auto* pts = ls->getCoordinatesRO();
+ if (!pts || pts->size() < 2) {
+ return;
+ }
+ DHDSegment prevSeg = DHDSegment::create(pts->getAt<CoordinateXY>(0),
+ pts->getAt<CoordinateXY>(1),
+ targetDistance);
+ consider(prevSeg);
+ for (std::size_t i = 1; i < pts->size() - 1; i++) {
+ DHDSegment seg = DHDSegment::create(prevSeg, pts->getAt<CoordinateXY>(i + 1), targetDistance);
+ consider(seg);
+ prevSeg = seg;
+ }
+ }
+
+ void addNonInterior(const DHDSegment& segment) const
+ {
+ if (segment.getMaxDistance() > 0.0) {
+ segQueue.push(segment);
+ return;
+ }
+ if (targetDistance.isInterior(segment.getEndpoint(0), segment.getEndpoint(1))) {
+ return;
+ }
+ segQueue.push(segment);
+ }
+
+ void consider(const DHDSegment& seg)
+ {
+ if (!segMaxDist || seg.getMaxDistanceBound() > segMaxDist->getMaxDistance()) {
+ addNonInterior(seg);
+ }
+ if (!segMaxDist || seg.getMaxDistance() > segMaxDist->getMaxDistance()) {
+ segMaxDist = seg;
+ }
+ }
+
+ std::priority_queue<DHDSegment>& segQueue;
+ TargetDistance& targetDistance;
+ std::optional<DHDSegment> segMaxDist;
+ };
+
+ LineFilter filter(segQueue, *targetDistance);
+ geom.apply_ro(&filter);
+
+ std::optional<DHDSegment> segMaxDist;
+ while (!segQueue.empty()) {
+ DHDSegment segMaxBound = segQueue.top();
+ segQueue.pop();
+
+ // Save if segment point is farther than current farthest
+ if (!segMaxDist || segMaxBound.getMaxDistance() > segMaxDist->getMaxDistance()) {
+ segMaxDist = segMaxBound;
+ }
+
+ // 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, isFullyWithin 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
+ && targetDistance->isSameOrCollinear(
+ segMaxBound.getEndpoint(0), segMaxBound.getEndpoint(1))) {
+ continue;
+ }
+
+ // 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) {
+ auto bisects = segMaxBound.bisect(*targetDistance);
+ filter.addNonInterior(bisects[0]);
+ filter.addNonInterior(bisects[1]);
+ }
+ }
+
+ // A segment at maximum distance was found.
+ // Return the farthest point pair
+ if (segMaxDist) {
+ 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
+ const CoordinateXY *maxPt = geom.getCoordinate();
+ if (!maxPt) {
+ return std::nullopt;
+ }
+ return pair(*maxPt, *maxPt);
+}
+
+std::optional<DirectedHausdorffDistance::PointPair>
+DirectedHausdorffDistance::computeForAreaInterior(
+ const Geometry& geom, double tolerance) const
+{
+ if (tolerance <= 0.0) {
+ return std::nullopt;
+ }
+ const Geometry& polygonal = geom;
+ if (polygonal.getEnvelopeInternal()->disjoint(target.getEnvelopeInternal())) {
+ return std::nullopt;
+ }
+
+ LargestEmptyCircle lec(&target, &polygonal, tolerance * AREA_INTERIOR_TOLERANCE_FACTOR);
+ auto centerPt = lec.getCenter();
+ const CoordinateXY* ptA = centerPt->getCoordinate();
+ if (!ptA) {
+ return std::nullopt;
+ }
+ if (targetDistance->isInterior(*ptA)) {
+ return std::nullopt;
+ }
+ CoordinateXY ptB = targetDistance->nearestFacetPoint(*ptA);
+ return pair(*ptA, ptB);
+}
+
+} // namespace distance
+} // namespace algorithm
+} // namespace geos
diff --git a/src/operation/distance/CoordinateSequenceLocation.cpp b/src/operation/distance/CoordinateSequenceLocation.cpp
new file mode 100644
index 000000000..7e6281a02
--- /dev/null
+++ b/src/operation/distance/CoordinateSequenceLocation.cpp
@@ -0,0 +1,107 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2026 Martin Davis
+ *
+ * This is free software; you can redistribute and/or modify it under
+ * the terms of the GNU Lesser General Public Licence as published
+ * by the Free Software Foundation.
+ * See the COPYING file for more information.
+ *
+ **********************************************************************
+ *
+ * Last port: operation/distance/CoordinateSequenceLocation.java (aff11591)
+ *
+ **********************************************************************/
+
+#include <geos/geom/CoordinateSequence.h>
+#include <geos/operation/distance/CoordinateSequenceLocation.h>
+
+using geos::geom::CoordinateXY;
+using geos::geom::CoordinateSequence;
+
+namespace geos {
+namespace operation {
+namespace distance {
+
+CoordinateSequenceLocation::CoordinateSequenceLocation(
+ const CoordinateSequence* p_seq, std::size_t p_index, const CoordinateXY& p_pt)
+ : seq(p_seq)
+ , index(p_index)
+ , pt(p_pt)
+{
+ if (seq && index >= seq->size() && seq->size() > 0) {
+ index = seq->size() - 1;
+ }
+}
+
+const CoordinateXY&
+CoordinateSequenceLocation::getCoordinate() const
+{
+ return pt;
+}
+
+std::size_t
+CoordinateSequenceLocation::getIndex() const
+{
+ return index;
+}
+
+bool
+CoordinateSequenceLocation::isSameSegment(const CoordinateSequenceLocation& f) const
+{
+ 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)) {
+ const CoordinateXY& endPt = seq->getAt<CoordinateXY>(index + 1);
+ return f.pt.equals2D(endPt);
+ }
+ if (isNext(f.index, index)) {
+ const CoordinateXY& endPt = f.seq->getAt<CoordinateXY>(index + 1);
+ return pt.equals2D(endPt);
+ }
+ return false;
+}
+
+bool
+CoordinateSequenceLocation::isNext(std::size_t p_index, std::size_t index1) const
+{
+ if (index1 == p_index + 1) {
+ return true;
+ }
+ // JTS writes `index1 == 0 && isRing && index1 == size-1`, which is
+ // unreachable. The ring wrap is last segment (size-2) adjacent to 0.
+ if (seq && seq->isRing() && index1 == 0 && p_index + 1 == seq->size() - 1) {
+ return true;
+ }
+ return false;
+}
+
+const CoordinateXY&
+CoordinateSequenceLocation::getEndPoint(int i) const
+{
+ if (i == 0) {
+ return seq->getAt<CoordinateXY>(index);
+ }
+ return seq->getAt<CoordinateXY>(index + 1);
+}
+
+std::size_t
+CoordinateSequenceLocation::normalize(std::size_t p_index) const
+{
+ if (seq && p_index >= seq->size() - 1 && seq->isRing()) {
+ return 0;
+ }
+ return p_index;
+}
+
+} // namespace distance
+} // namespace operation
+} // namespace geos
diff --git a/src/operation/distance/FacetSequence.cpp b/src/operation/distance/FacetSequence.cpp
index ebc1501dd..2b08f1b72 100644
--- a/src/operation/distance/FacetSequence.cpp
+++ b/src/operation/distance/FacetSequence.cpp
@@ -12,10 +12,11 @@
*
**********************************************************************
*
- * Last port: operation/distance/FacetSequence.java (f6187ee2 JTS-1.14)
+ * Last port: operation/distance/FacetSequence.java (aff11591)
*
**********************************************************************/
+#include <geos/geom/CoordinateSequence.h>
#include <geos/geom/Geometry.h>
#include <geos/geom/LineSegment.h>
#include <geos/algorithm/Distance.h>
@@ -229,3 +230,57 @@ FacetSequence::getCoordinate(std::size_t index) const
return &(pts->getAt(start + index));
}
+std::size_t
+FacetSequence::normalize(const CoordinateSequence& pts, std::size_t index)
+{
+ if (index >= pts.size() - 1 && pts.isRing()) {
+ return 0;
+ }
+ return index;
+}
+
+CoordinateSequenceLocation
+FacetSequence::nearestLocation(const CoordinateXY& p) const
+{
+ if (isPoint()) {
+ // JTS uses index 0 / pts[0]; use start so a one-point facet
+ // that is a subsequence of a longer CoordinateSequence is correct.
+ return CoordinateSequenceLocation(pts, start, pts->getAt(start));
+ }
+ return nearestLocationOnLine(p);
+}
+
+CoordinateSequenceLocation
+FacetSequence::nearestLocationOnLine(const CoordinateXY& pt) const
+{
+ double minDistance = DoubleInfinity;
+ std::size_t index = start;
+ Coordinate nearestPt;
+
+ const Coordinate queryPt(pt);
+ for (std::size_t i = start; i < end - 1; i++) {
+ const Coordinate& q0 = pts->getAt(i);
+ const Coordinate& q1 = pts->getAt(i + 1);
+ double dist = Distance::pointToSegment(queryPt, q0, q1);
+ if (dist < minDistance) {
+ minDistance = dist;
+ LineSegment seg(q0, q1);
+ seg.closestPoint(queryPt, nearestPt);
+ 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 && queryPt.equals2D(q1)) {
+ if (index < pts->size() - 1) {
+ index++;
+ }
+ //-- normalize index for a ring
+ index = normalize(*pts, index);
+ }
+ if (minDistance <= 0.0) {
+ break;
+ }
+ }
+ }
+ return CoordinateSequenceLocation(pts, index, nearestPt);
+}
+
diff --git a/src/operation/distance/IndexedFacetDistance.cpp b/src/operation/distance/IndexedFacetDistance.cpp
index 6555920cd..babe15438 100644
--- a/src/operation/distance/IndexedFacetDistance.cpp
+++ b/src/operation/distance/IndexedFacetDistance.cpp
@@ -12,13 +12,16 @@
*
**********************************************************************
*
- * Last port: operation/distance/IndexedFacetDistance.java (f6187ee2 JTS-1.14)
+ * Last port: operation/distance/IndexedFacetDistance.java (aff11591)
*
**********************************************************************/
#include <geos/geom/Coordinate.h>
+#include <geos/geom/CoordinateSequence.h>
#include <geos/index/strtree/STRtree.h>
#include <geos/operation/distance/IndexedFacetDistance.h>
+#include <geos/operation/distance/FacetSequence.h>
+#include <geos/util/GEOSException.h>
using namespace geos::geom;
using namespace geos::index::strtree;
@@ -105,6 +108,48 @@ IndexedFacetDistance::isWithinDistance(const Geometry* g, double maxDistance) co
return cachedTree->isWithinDistance<FacetDistance>(*tree2, maxDistance);
}
+CoordinateSequenceLocation
+IndexedFacetDistance::nearestLocation(const geom::CoordinateXY& p) const
+{
+ // JTS wraps p in a CoordinateArraySequence; GEOS CoordinateSequence is the equivalent.
+ CoordinateSequence seq{CoordinateXY(p.x, p.y)};
+ FacetSequence query(&seq, 0, 1);
+ FacetDistance itemDist;
+ const FacetSequence* nearest = cachedTree->nearestNeighbour<FacetDistance>(
+ *query.getEnvelope(), &query, itemDist);
+ if (!nearest) {
+ throw util::GEOSException("Cannot calculate IndexedFacetDistance on empty geometries.");
+ }
+ return nearest->nearestLocation(p);
+}
+
+geom::CoordinateXY
+IndexedFacetDistance::nearestPoint(const geom::CoordinateXY& p) const
+{
+ return nearestLocation(p).getCoordinate();
+}
+
+double
+IndexedFacetDistance::distance(const geom::CoordinateXY& p) const
+{
+ return p.distance(nearestPoint(p));
+}
+
+double
+IndexedFacetDistance::distance(const geom::CoordinateXY& p0, const geom::CoordinateXY& p1) const
+{
+ CoordinateSequence seq{CoordinateXY(p0.x, p0.y), CoordinateXY(p1.x, p1.y)};
+ FacetSequence query(&seq, 0, 2);
+ FacetDistance itemDist;
+ const FacetSequence* nearest = cachedTree->nearestNeighbour<FacetDistance>(
+ *query.getEnvelope(), &query, itemDist);
+ if (!nearest) {
+ throw util::GEOSException("Cannot calculate IndexedFacetDistance on empty geometries.");
+ }
+ auto locs = nearest->nearestLocations(query);
+ return locs[0].distance(locs[1]);
+}
+
}
}
diff --git a/tests/unit/algorithm/distance/DirectedHausdorffDistanceTest.cpp b/tests/unit/algorithm/distance/DirectedHausdorffDistanceTest.cpp
new file mode 100644
index 000000000..5d9c4ce09
--- /dev/null
+++ b/tests/unit/algorithm/distance/DirectedHausdorffDistanceTest.cpp
@@ -0,0 +1,630 @@
+//
+// Test Suite for geos::algorithm::distance::DirectedHausdorffDistance
+// Ported from JTS DirectedHausdorffDistanceTest (locationtech/jts#1182)
+// plus the discrete-under-estimate witness from DiscreteHausdorffDistance.
+
+#include <tut/tut.hpp>
+#include <tut/tut_macros.hpp>
+
+#include <geos/algorithm/distance/DirectedHausdorffDistance.h>
+#include <geos/algorithm/distance/DiscreteHausdorffDistance.h>
+#include <geos/geom/Coordinate.h>
+#include <geos/geom/CoordinateSequence.h>
+#include <geos/geom/Geometry.h>
+#include <geos/geom/GeometryFactory.h>
+#include <geos/geom/PrecisionModel.h>
+#include <geos/io/WKTReader.h>
+#include <geos/util/IllegalArgumentException.h>
+
+#include <cmath>
+#include <memory>
+#include <string>
+
+#include "utility.h"
+
+using geos::algorithm::distance::DirectedHausdorffDistance;
+using geos::algorithm::distance::DiscreteHausdorffDistance;
+using geos::geom::CoordinateXY;
+using geos::geom::Geometry;
+using geos::geom::GeometryFactory;
+using geos::geom::PrecisionModel;
+
+namespace tut {
+
+struct test_directedhausdorffdistance_data {
+ test_directedhausdorffdistance_data()
+ : pm()
+ , gf(GeometryFactory::create(&pm))
+ , reader(gf.get())
+ {}
+
+ static constexpr double TOLERANCE = 0.001;
+
+ std::unique_ptr<Geometry>
+ read(const std::string& wkt) const
+ {
+ return reader.read(wkt);
+ }
+
+ void
+ checkDistance(const std::string& wkt1, const std::string& wkt2,
+ double expectedDistance) const
+ {
+ auto g1 = read(wkt1);
+ auto g2 = read(wkt2);
+
+ double dist = DirectedHausdorffDistance::distance(*g1, *g2);
+ ensure(std::fabs(dist - expectedDistance) <= TOLERANCE);
+ }
+
+ void checkDistance(const std::string& wkt1, const std::string& wkt2,
+ double tolerance, const std::string& wktExpected) const
+ {
+ auto g1 = read(wkt1);
+ auto g2 = read(wkt2);
+
+ auto pts = DirectedHausdorffDistance::distancePoints(*g1, *g2, tolerance);
+ ensure(pts.has_value());
+
+ auto seq = std::make_unique<geos::geom::CoordinateSequence>(2);
+ seq->setAt(pts.value()[0], 0);
+ seq->setAt(pts.value()[1], 1);
+
+ auto result = g1->getFactory()->createLineString(std::move(seq));
+ ensure_equals_exact_geometry(static_cast<const Geometry*>(result.get()), read(wktExpected).get(), TOLERANCE);
+ }
+
+ void checkDistanceStartPtLen(const std::string wkt1, const std::string& wkt2,
+ const std::string& wktExpected, double resultTolerance) const
+ {
+ auto g1 = read(wkt1);
+ auto g2 = read(wkt2);
+
+ auto pts = DirectedHausdorffDistance::distancePoints(*g1, *g2);
+ ensure(pts.has_value());
+
+ auto seq = std::make_unique<geos::geom::CoordinateSequence>(2);
+ seq->setAt(pts.value()[0], 0);
+ seq->setAt(pts.value()[1], 1);
+
+ auto result = g1->getFactory()->createLineString(std::move(seq));
+ auto expected = reader.read<LineString>(wktExpected);
+
+ const CoordinateXY& resultPt = result->getCoordinatesRO()->getAt<CoordinateXY>(0);
+ const CoordinateXY& expectedPt = expected->getCoordinatesRO()->getAt<CoordinateXY>(0);
+
+ ensure_equals_xy(expectedPt, resultPt, resultTolerance);
+
+ double distResult = result->getLength();
+ double distExpected = expected->getLength();
+ ensure_equals("distance is not within tolerance of expected", distExpected, distResult, resultTolerance);
+ }
+
+ void
+ checkDistance(const std::string& wkt1, const std::string& wkt2,
+ double tolerance, double expectedDistance) const
+ {
+ auto g1 = read(wkt1);
+ auto g2 = read(wkt2);
+ double dist = DirectedHausdorffDistance::distance(*g1, *g2, tolerance);
+ ensure(std::fabs(dist - expectedDistance) <= TOLERANCE);
+ }
+
+ void
+ checkDistance(const std::string& wkt1, const std::string& wkt2,
+ const std::string& wktExpected) const
+ {
+ auto g1 = read(wkt1);
+ auto g2 = read(wkt2);
+ auto pts = DirectedHausdorffDistance::distancePoints(*g1, *g2);
+ ensure(pts.has_value());
+
+ auto seq = std::make_unique<geos::geom::CoordinateSequence>(2);
+ seq->setAt(pts.value()[0], 0);
+ seq->setAt(pts.value()[1], 1);
+
+ auto result = g1->getFactory()->createLineString(std::move(seq));
+ ensure_equals_exact_geometry(static_cast<const Geometry*>(result.get()), read(wktExpected).get(), TOLERANCE);
+
+ ensure_equals_exact_geometry(static_cast<const Geometry*>(result.get()), read(wktExpected).get(), TOLERANCE);
+ }
+
+ void
+ checkHausdorff(const std::string& wkt1, const std::string& wkt2,
+ const std::string& wktExpected) const
+ {
+ std::unique_ptr<Geometry> g1 = read(wkt1);
+ std::unique_ptr<Geometry> g2 = read(wkt2);
+
+ auto pts = DirectedHausdorffDistance::hausdorffDistancePoints(*g1, *g2);
+ ensure(pts.has_value());
+
+ auto seq = std::make_unique<geos::geom::CoordinateSequence>(2);
+ seq->setAt(pts.value()[0], 0);
+ seq->setAt(pts.value()[1], 1);
+
+ auto result = g1->getFactory()->createLineString(std::move(seq));
+ ensure_equals_exact_geometry(static_cast<const Geometry*>(result.get()), read(wktExpected).get(), TOLERANCE);
+ }
+
+ void
+ checkDistanceEmpty(const std::string& a, const std::string& b) const
+ {
+ std::unique_ptr<Geometry> g1 = read(a);
+ std::unique_ptr<Geometry> g2 = read(b);
+ ensure(!DirectedHausdorffDistance::distancePoints(*g1, *g2).has_value());
+ ensure(std::isnan(DirectedHausdorffDistance::distance(*g1, *g2)));
+ ensure(std::isnan(DirectedHausdorffDistance::hausdorffDistance(*g1, *g2)));
+ }
+
+ void
+ checkFullyWithinDistance(const std::string& a, const std::string& b,
+ double distance, bool expected) const
+ {
+ std::unique_ptr<Geometry> g1 = read(a);
+ std::unique_ptr<Geometry> g2 = read(b);
+ bool result = DirectedHausdorffDistance::isFullyWithinDistance(*g1, *g2, distance);
+ ensure_equals(result, expected);
+ }
+
+ void checkFullyWithinDistanceEmpty(const std::string& a, const std::string& b) const
+ {
+ 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);
+ }
+
+ PrecisionModel pm;
+ GeometryFactory::Ptr gf;
+ geos::io::WKTReader reader;
+};
+
+typedef test_group<test_directedhausdorffdistance_data> group;
+typedef group::object object;
+
+group test_DirectedHausdorffDistance_group(
+ "geos::algorithm::distance::DirectedHausdorffDistance");
+
+template<>
+template<>
+void object::test<1>()
+{
+ set_test_name("testEmpty");
+
+ checkDistanceEmpty("POINT EMPTY", "POINT (1 1)");
+ checkDistanceEmpty("LINESTRING EMPTY", "LINESTRING (0 0, 2 1)");
+ checkDistanceEmpty("POLYGON EMPTY", "POLYGON ((1 9, 9 9, 9 1, 1 1, 1 9))");
+}
+
+template<>
+template<>
+void object::test<2>()
+{
+ set_test_name("testZeroTolerancePoint");
+
+ checkDistance("POINT (5 5)", "LINESTRING (5 1, 9 5)",
+ 0,
+ "LINESTRING (5 5, 7 3)");
+}
+
+template<>
+template<>
+void object::test<3>()
+{
+ set_test_name("testZeroToleranceLine");
+
+ checkDistance("LINESTRING (1 5, 5 5)", "LINESTRING (5 1, 9 5)",
+ 0,
+ "LINESTRING (1 5, 5 1)");
+}
+
+template<>
+template<>
+void object::test<4>()
+{
+ set_test_name("testZeroToleranceZeroLengthLineQuery");
+
+ checkDistance("LINESTRING (5 5, 5 5)", "LINESTRING (5 1, 9 5)",
+ 0,
+ "LINESTRING (5 5, 7 3)");
+}
+
+template<>
+template<>
+void object::test<5>()
+{
+ set_test_name("testZeroLengthLineQuery");
+
+ checkDistance("LINESTRING (5 5, 5 5)", "LINESTRING (5 1, 9 5)",
+ "LINESTRING (5 5, 7 3)");
+}
+
+template<>
+template<>
+void object::test<6>()
+{
+ set_test_name("testZeroLengthPolygonQuery");
+
+ checkDistance("POLYGON ((5 5, 5 5, 5 5, 5 5))", "LINESTRING (5 1, 9 5)",
+ "LINESTRING (5 5, 7 3)");
+}
+
+template<>
+template<>
+void object::test<7>()
+{
+ set_test_name("testZeroLengthLineTarget");
+
+ checkDistance("POINT (5 5)", "LINESTRING (5 1, 5 1)",
+ "LINESTRING (5 5, 5 1)");
+}
+
+template<>
+template<>
+void object::test<8>()
+{
+ set_test_name("testNegativeTolerancePoint");
+
+ auto g1 = read("POINT (5 5)");
+ auto g2 = read("LINESTRING (5 1, 9 5)");
+
+ ensure_THROW(DirectedHausdorffDistance::distance(*g1, *g2, -1.0), geos::util::IllegalArgumentException);
+}
+
+template<>
+template<>
+void object::test<9>()
+{
+ set_test_name("testNegativeToleranceLine");
+
+ auto g1 = read("LINESTRING (1 1, 5 5)");
+ auto g2 = read("LINESTRING (5 1, 9 5)");
+
+ ensure_THROW(DirectedHausdorffDistance::distance(*g1, *g2, -1.0), geos::util::IllegalArgumentException);
+}
+
+template<>
+template<>
+void object::test<10>()
+{
+ set_test_name("testPointPoint");
+
+ checkHausdorff("POINT (0 0)", "POINT (1 1)", "LINESTRING (0 0, 1 1)");
+}
+
+template<>
+template<>
+void object::test<11>()
+{
+ set_test_name("testPointPoints");
+
+ std::string a = "MULTIPOINT ((0 1), (2 3), (4 5), (6 6))";
+ std::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)");
+}
+
+template<>
+template<>
+void object::test<12>()
+{
+ set_test_name("testPointPolygonInterior");
+
+ checkDistance("POINT (3 4)", "POLYGON ((1 9, 9 9, 9 1, 1 1, 1 9))",
+ 0);
+}
+
+template<>
+template<>
+void object::test<13>()
+{
+ set_test_name("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)");
+}
+
+template<>
+template<>
+void object::test<14>()
+{
+ set_test_name("testLineSegments");
+
+ checkHausdorff("LINESTRING (0 0, 2 0)", "LINESTRING (0 0, 2 1)",
+ "LINESTRING (2 0, 2 1)");
+}
+
+template<>
+template<>
+void object::test<15>()
+{
+ set_test_name("testLineSegments2");
+
+ checkHausdorff("LINESTRING (0 0, 2 0)", "LINESTRING (0 1, 1 2, 2 1)",
+ "LINESTRING (1 0, 1 2)");
+}
+
+template<>
+template<>
+void object::test<16>()
+{
+ set_test_name("testLinePoints");
+
+ checkHausdorff("LINESTRING (0 0, 2 0)", "MULTIPOINT (0 2, 1 0, 2 1)",
+ "LINESTRING (0 0, 0 2)");
+}
+
+template<>
+template<>
+void object::test<17>()
+{
+ set_test_name("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);
+}
+
+template<>
+template<>
+void object::test<18>()
+{
+ set_test_name("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)");
+}
+
+template<>
+template<>
+void object::test<19>()
+{
+ set_test_name("testLinesPolygon2");
+
+ std::string a = "MULTILINESTRING ((2 3, 2 7), (9 1, 9 8, 4 9))";
+ std::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)");
+}
+
+template<>
+template<>
+void object::test<20>()
+{
+ set_test_name("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)");
+}
+
+template<>
+template<>
+void object::test<21>()
+{
+ set_test_name("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);
+}
+
+template<>
+template<>
+void object::test<22>()
+{
+ set_test_name("testPolygonPolygon");
+
+ std::string a = "POLYGON ((2 18, 18 18, 17 3, 2 2, 2 18))";
+ std::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)");
+}
+
+template<>
+template<>
+void object::test<23>()
+{
+ set_test_name("testPolygonPolygonHolesNested");
+
+ // B is contained in A
+ std::string a = "POLYGON ((1 19, 19 19, 19 1, 1 1, 1 19), (6 8, 11 14, 15 7, 6 8))";
+ std::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);
+}
+
+template<>
+template<>
+void object::test<24>()
+{
+ set_test_name("testMultiPolygons");
+
+ std::string a = "MULTIPOLYGON (((1 1, 1 10, 5 1, 1 1)), ((4 17, 9 15, 9 6, 4 17)))";
+ std::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);
+}
+
+template<>
+template<>
+void object::test<25>()
+{
+ set_test_name("testLinePolygonCrossing");
+
+ std::string wkt1 = "LINESTRING (2 5, 5 10, 6 4)";
+ std::string wkt2 = "POLYGON ((1 9, 9 9, 9 1, 1 1, 1 9))";
+ checkDistance(wkt1, wkt2, "LINESTRING (5 10, 5 9)");
+}
+
+template<>
+template<>
+void object::test<26>()
+{
+ set_test_name("testNonVertexResult");
+
+ std::string wkt1 = "LINESTRING (1 1, 5 10, 9 1)";
+ std::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)");
+}
+
+template<>
+template<>
+void object::test<27>()
+{
+ set_test_name("testDirectedLines");
+
+ std::string wkt1 = "LINESTRING (1 6, 3 5, 1 4)";
+ std::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)");
+}
+
+template<>
+template<>
+void object::test<28>()
+{
+ set_test_name("testDirectedLines2");
+
+ std::string wkt1 = "LINESTRING (1 6, 3 5, 1 4)";
+ std::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.
+ */
+template<>
+template<>
+void object::test<29>()
+{
+ set_test_name("testInteriorSegmentsLargeTol");
+
+ std::string a = "POLYGON ((4 6, 5 6, 5 5, 4 5, 4 6))";
+ std::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
+ */
+template<>
+template<>
+void object::test<30>()
+{
+ set_test_name("testInteriorSegmentsSameExterior");
+ std::string a = "POLYGON ((1 9, 3 9, 4 5, 5.05 9, 9 9, 9 1, 1 1, 1 9))";
+ std::string b = "POLYGON ((1 9, 9 9, 9 1, 1 1, 1 9))";
+ checkDistance(a, b, 0.0);
+}
+
+//-----------------------------------------------------
+
+template<>
+template<>
+void object::test<31>()
+{
+ set_test_name("testFullyWithinDistanceEmptyPoints");
+ std::string a = "POINT EMPTY";
+ std::string b = "MULTIPOINT ((1 1), (9 9))";
+ checkFullyWithinDistanceEmpty(a, b);
+}
+
+template<>
+template<>
+void object::test<32>()
+{
+ set_test_name("testFullyWithinDistanceEmptyLine");
+ std::string a = "LINESTRING EMPTY";
+ std::string b = "LINESTRING (9 9, 1 1)";
+ checkFullyWithinDistanceEmpty(a, b);
+}
+
+template<>
+template<>
+void object::test<33>()
+{
+ //-- shows withinDistance envelope check not triggering for disconnected A
+ set_test_name("testFullyWithinDistancePoints");
+ std::string a = "MULTIPOINT ((1 9), (9 1))";
+ std::string b = "MULTIPOINT ((1 1), (9 9))";
+ checkFullyWithinDistance(a, b, 1, false);
+ checkFullyWithinDistance(a, b, 8.1, true);
+}
+
+template<>
+template<>
+void object::test<34>()
+{
+ set_test_name("testFullyWithinDistanceDisconnectedLines");
+ std::string a = "MULTILINESTRING ((1 9, 2 9), (8 1, 9 1))";
+ std::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);
+}
+
+template<>
+template<>
+void object::test<35>()
+{
+ set_test_name("testFullyWithinDistanceDisconnectedPolygons");
+ std::string a = "MULTIPOLYGON (((1 9, 2 9, 2 8, 1 8, 1 9)), ((8 2, 9 2, 9 1, 8 1, 8 2)))";
+ std::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);
+}
+
+template<>
+template<>
+void object::test<36>()
+{
+ set_test_name("testFullyWithinDistanceLines");
+ std::string a = "MULTILINESTRING ((1 1, 3 3), (7 7, 9 9))";
+ std::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);
+}
+
+template<>
+template<>
+void object::test<37>()
+{
+ set_test_name("testFullyWithinDistancePolygons");
+ std::string a = "POLYGON ((1 4, 4 4, 4 1, 1 1, 1 4))";
+ std::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);
+}
+
+template<>
+template<>
+void object::test<38>()
+{
+ set_test_name("testFullyWithinDistancePolygonsNestedWithHole");
+ std::string a = "POLYGON ((2 8, 8 8, 8 2, 2 2, 2 8))";
+ std::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);
+}
+
+
+} // namespace tut
diff --git a/tests/unit/capi/GEOSDirectedHausdorffDistanceTest.cpp b/tests/unit/capi/GEOSDirectedHausdorffDistanceTest.cpp
new file mode 100644
index 000000000..3832542e9
--- /dev/null
+++ b/tests/unit/capi/GEOSDirectedHausdorffDistanceTest.cpp
@@ -0,0 +1,96 @@
+//
+// Test Suite for C-API GEOSDirectedHausdorffDistance /
+// GEOSSymmetricHausdorffDistance
+
+#include <tut/tut.hpp>
+#include <geos_c.h>
+
+#include "capi_test_utils.h"
+
+#include <cmath>
+
+namespace tut {
+
+struct test_capigeosdirectedhausdorffdistance_data : public capitest::utility {
+};
+
+typedef test_group<test_capigeosdirectedhausdorffdistance_data> group;
+typedef group::object object;
+
+group test_capigeosdirectedhausdorffdistance_group(
+ "capi::GEOSDirectedHausdorffDistance");
+
+template<>
+template<>
+void object::test<1>()
+{
+ set_test_name("maximum distance point pair within segment interior");
+
+ geom1_ = fromWKT("LINESTRING (0 0, 100 0, 10 100, 10 100)");
+ geom2_ = fromWKT("LINESTRING (0 100, 0 10, 80 10)");
+
+ double discrete = 0.0;
+ double directed12 = 0.0;
+ double directed21 = 0.0;
+ double symmetric = 0.0;
+ ensure_equals(GEOSHausdorffDistance(geom1_, geom2_, &discrete), 1);
+ ensure_equals(GEOSDirectedHausdorffDistance(geom1_, geom2_, &directed12), 1);
+ ensure_equals(GEOSDirectedHausdorffDistance(geom2_, geom1_, &directed21), 1);
+ ensure_equals(GEOSSymmetricHausdorffDistance(geom1_, geom2_, &symmetric), 1);
+
+ ensure_equals("GEOSHausdorffDistance", discrete, 22.36, 1e-2);
+ ensure_equals("Directed 1->2", directed12, 47.891845703125);
+ ensure_equals("Directed 2->1", directed21, 44.533046060947164);
+ ensure_equals("Symmetric", symmetric, 47.891845703125);
+}
+
+template<>
+template<>
+void object::test<2>()
+{
+ set_test_name("empty writes NaN, not 0");
+
+ geom1_ = fromWKT("LINESTRING EMPTY");
+ geom2_ = fromWKT("LINESTRING (0 0, 2 1)");
+
+ double dist = 0.0;
+ ensure_equals(GEOSDirectedHausdorffDistance(geom1_, geom2_, &dist), 1);
+ ensure(std::isnan(dist));
+}
+
+template<>
+template<>
+void object::test<3>()
+{
+ set_test_name("GEOSDirectedHausdorffDistanceWithPoints");
+
+ geom1_ = fromWKT("POINT (0 0)");
+ geom2_ = fromWKT("POINT (3 4)");
+
+ double dist = 0.0;
+ double p1x = 0, p1y = 0, p2x = 0, p2y = 0;
+ ensure_equals(GEOSDirectedHausdorffDistanceWithPoints(
+ geom1_, geom2_, &dist, &p1x, &p1y, &p2x, &p2y), 1);
+ ensure_distance(dist, 5.0, 1e-12);
+ ensure_equals(p1x, 0.0);
+ ensure_equals(p1y, 0.0);
+ ensure_equals(p2x, 3.0);
+ ensure_equals(p2y, 4.0);
+}
+
+template<>
+template<>
+void object::test<4>()
+{
+ set_test_name("GEOSDirectedHausdorffDistanceWithin");
+
+ geom1_ = fromWKT("LINESTRING (0 0, 2 0)");
+ geom2_ = fromWKT("LINESTRING (0 0, 2 1)");
+
+ ensure_equals(GEOSDirectedHausdorffDistanceWithin(geom1_, geom2_, 0.5), 0);
+ ensure_equals(GEOSDirectedHausdorffDistanceWithin(geom1_, geom2_, 1.0), 1);
+ ensure_equals(GEOSDirectedHausdorffDistanceWithin(geom1_, geom2_, 2.0), 1);
+}
+
+
+} // namespace tut
\ No newline at end of file
diff --git a/util/geosop/GeometryOp.cpp b/util/geosop/GeometryOp.cpp
index 29e0dd59f..0c0bcfb34 100644
--- a/util/geosop/GeometryOp.cpp
+++ b/util/geosop/GeometryOp.cpp
@@ -32,6 +32,7 @@
#include <geos/algorithm/MinimumDiameter.h>
#include <geos/algorithm/MinimumBoundingCircle.h>
#include <geos/algorithm/distance/DiscreteHausdorffDistance.h>
+#include <geos/algorithm/distance/DirectedHausdorffDistance.h>
#include <geos/algorithm/distance/DiscreteFrechetDistance.h>
#include <geos/algorithm/hull/ConcaveHull.h>
#include <geos/geom/util/Densifier.h>
@@ -643,6 +644,22 @@ std::vector<GeometryOpCreator> opRegistry {
return new Result( geos::algorithm::distance::DiscreteHausdorffDistance::distance(geom, geomB ) );
});
}},
+{"directedHausdorffDistance", [](std::string name) { return GeometryOp::create(name,
+ catDist,
+ "compute directed (locus) Hausdorff distance from geometry A to B",
+ Result::typeDouble,
+ [](const Geometry& geom, const Geometry& geomB) {
+ return new Result( geos::algorithm::distance::DirectedHausdorffDistance::distance(geom, geomB ) );
+ });
+ }},
+{"symmetricHausdorffDistance", [](std::string name) { return GeometryOp::create(name,
+ catDist,
+ "compute symmetric (locus) Hausdorff distance between geometry A and B",
+ Result::typeDouble,
+ [](const Geometry& geom, const Geometry& geomB) {
+ return new Result( geos::algorithm::distance::DirectedHausdorffDistance::hausdorffDistance(geom, geomB ) );
+ });
+ }},
/*
// MD - can't get this to work for now
add("frechetDistanceLine", 2, 0, Result::typeGeometry, catDist,
-----------------------------------------------------------------------
Summary of changes:
NEWS.md | 4 +
Version.txt | 4 +-
capi/geos_c.cpp | 24 +
capi/geos_c.h.in | 116 ++++
capi/geos_ts_c.cpp | 63 ++
.../algorithm/distance/DirectedHausdorffDistance.h | 300 ++++++++++
.../distance/CoordinateSequenceLocation.h | 76 +++
include/geos/operation/distance/FacetSequence.h | 21 +-
.../geos/operation/distance/IndexedFacetDistance.h | 18 +-
.../distance/DirectedHausdorffDistance.cpp | 659 +++++++++++++++++++++
.../distance/CoordinateSequenceLocation.cpp | 107 ++++
src/operation/distance/FacetSequence.cpp | 57 +-
src/operation/distance/IndexedFacetDistance.cpp | 47 +-
.../distance/DirectedHausdorffDistanceTest.cpp | 630 ++++++++++++++++++++
.../capi/GEOSDirectedHausdorffDistanceTest.cpp | 96 +++
util/geosop/GeometryOp.cpp | 17 +
16 files changed, 2231 insertions(+), 8 deletions(-)
create mode 100644 include/geos/algorithm/distance/DirectedHausdorffDistance.h
create mode 100644 include/geos/operation/distance/CoordinateSequenceLocation.h
create mode 100644 src/algorithm/distance/DirectedHausdorffDistance.cpp
create mode 100644 src/operation/distance/CoordinateSequenceLocation.cpp
create mode 100644 tests/unit/algorithm/distance/DirectedHausdorffDistanceTest.cpp
create mode 100644 tests/unit/capi/GEOSDirectedHausdorffDistanceTest.cpp
hooks/post-receive
--
GEOS
More information about the geos-commits
mailing list