From dd8fa1cd445db73a0a5e92333a23f29436490068 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Wed, 17 Dec 2025 00:35:58 +0100 Subject: [PATCH 1/2] DistanceOp::distance(): make it return NaN instead of 0 if one of the geometry is empty Fixes https://github.com/OSGeo/gdal/issues/12978 --- capi/geos_c.h.in | 7 ++++--- src/operation/distance/DistanceOp.cpp | 5 +++-- .../operation/distance/DistanceOpTest.cpp | 19 +++++++++---------- tests/xmltester/XMLTester.cpp | 8 ++++++-- .../xmltester/tests/general/TestDistance.xml | 12 ++++++------ 5 files changed, 28 insertions(+), 23 deletions(-) diff --git a/capi/geos_c.h.in b/capi/geos_c.h.in index 942663671e..78bb3b401f 100644 --- a/capi/geos_c.h.in +++ b/capi/geos_c.h.in @@ -3387,7 +3387,7 @@ extern char GEOS_DLL GEOSisSimple(const GEOSGeometry* g); * In one step, calculate and return whether a geometry is simple * and one more more points at which the geometry self-intersects * at interior points. -* Caller has the responsibility to destroy 'location' with +* Caller has the responsibility to destroy 'location' with * GEOSGeom_destroy() * * \param g The geometry to test @@ -3631,7 +3631,8 @@ extern int GEOS_DLL GEOSGeomGetLength( * Calculate the distance between two geometries. * \param[in] g1 Input geometry * \param[in] g2 Input geometry -* \param[out] dist Pointer to be filled in with distance result +* \param[out] dist Pointer to be filled in with distance result. NaN is returned +* if one of the geometry is empty. * \return 1 on success, 0 on exception. * \since 2.2 */ @@ -5462,7 +5463,7 @@ extern char GEOS_DLL GEOSIntersects(const GEOSGeometry* g1, const GEOSGeometry* extern char GEOS_DLL GEOSCrosses(const GEOSGeometry* g1, const GEOSGeometry* g2); /** -* Tests if geometry g1 is completely within g2, +* Tests if geometry g1 is completely within g2, * but not wholly contained in the boundary of g2. * \param g1 Input geometry * \param g2 Input geometry diff --git a/src/operation/distance/DistanceOp.cpp b/src/operation/distance/DistanceOp.cpp index c6da7508fa..cf457fbe30 100644 --- a/src/operation/distance/DistanceOp.cpp +++ b/src/operation/distance/DistanceOp.cpp @@ -36,6 +36,7 @@ #include #include +#include #include #include @@ -97,7 +98,7 @@ DistanceOp::DistanceOp(const Geometry& g0, const Geometry& g1, double tdist) /** * Report the distance between the closest points on the input geometries. * - * @return the distance between the geometries + * @return the distance between the geometries, or NaN if one of them is empty. */ double DistanceOp::distance() @@ -111,7 +112,7 @@ DistanceOp::distance() throw IllegalArgumentException("null geometries are not supported"); } if(geom[0]->isEmpty() || geom[1]->isEmpty()) { - return 0.0; + return std::numeric_limits::quiet_NaN(); } if(geom[0]->getGeometryTypeId() == GEOS_POINT && geom[1]->getGeometryTypeId() == GEOS_POINT) { return static_cast(geom[0])->getCoordinate()->distance(*static_cast(geom[1])->getCoordinate()); diff --git a/tests/unit/operation/distance/DistanceOpTest.cpp b/tests/unit/operation/distance/DistanceOpTest.cpp index 32eee5bbff..cbd96f8060 100644 --- a/tests/unit/operation/distance/DistanceOpTest.cpp +++ b/tests/unit/operation/distance/DistanceOpTest.cpp @@ -228,7 +228,7 @@ void object::test<8> DistanceOp dist(g0.get(), g1.get()); - ensure_equals(dist.distance(), 0); + ensure(std::isnan(dist.distance())); ensure(dist.nearestPoints() == nullptr); } @@ -416,7 +416,7 @@ void object::test<16> DistanceOp dist(g0.get(), g1.get()); - ensure_equals(dist.distance(), 0); + ensure(std::isnan(dist.distance())); ensure(dist.nearestPoints() == nullptr); } @@ -497,7 +497,7 @@ void object::test<19> ensure(g1->isValid()); ensure(g2->isValid()); - ensure_equals(g1->distance(g2.get()), 0); + ensure(std::isnan(g1->distance(g2.get()))); } // Test case reported in Shapely @@ -580,7 +580,6 @@ void object::test<22>() ensure_equals(g2->distance(g1.get()), 1); } -// Empty is same as empty so zero...? template<> template<> void object::test<23>() @@ -589,8 +588,8 @@ void object::test<23>() auto g2 = wktreader.read("LINESTRING EMPTY"); ensure(g1 != nullptr && g2 != nullptr); - ensure_equals(g1->distance(g2.get()), 0); - ensure_equals(g2->distance(g1.get()), 0); + ensure(std::isnan(g1->distance(g2.get()))); + ensure(std::isnan(g2->distance(g1.get()))); } template<> @@ -601,8 +600,8 @@ void object::test<24>() auto g2 = wktreader.read("LINESTRING EMPTY"); ensure(g1 != nullptr && g2 != nullptr); - ensure_equals(g1->distance(g2.get()), 0); - ensure_equals(g2->distance(g1.get()), 0); + ensure(std::isnan(g1->distance(g2.get()))); + ensure(std::isnan(g2->distance(g1.get()))); } // But ignore empty if there's a real distance? @@ -638,8 +637,8 @@ void object::test<27>() auto g2 = wktreader.read("GEOMETRYCOLLECTION(POINT(1 0))"); ensure(g1 != nullptr && g2 != nullptr); - ensure_equals(g1->distance(g2.get()), 0); - ensure_equals(g2->distance(g1.get()), 0); + ensure(std::isnan(g1->distance(g2.get()))); + ensure(std::isnan(g2->distance(g1.get()))); } diff --git a/tests/xmltester/XMLTester.cpp b/tests/xmltester/XMLTester.cpp index c3b33b9394..acdba782da 100644 --- a/tests/xmltester/XMLTester.cpp +++ b/tests/xmltester/XMLTester.cpp @@ -894,8 +894,9 @@ Test::checkResult( bool result ) void Test::checkResult( double result) { - char* rest; - double expectedRes = std::strtod(opResult.c_str(), &rest); + char* rest = nullptr; + const double expectedRes = (opResult == "NaN" || opResult == "nan") ? + std::numeric_limits::quiet_NaN() : std::strtod(opResult.c_str(), &rest); if(rest == opResult.c_str()) { throw std::runtime_error("malformed testcase: missing expected double value"); } @@ -904,6 +905,9 @@ Test::checkResult( double result) isSuccess = true; } } + else if (std::isnan(expectedRes)) { + isSuccess = std::isnan(result); + } else { if (std::abs(expectedRes - result) / expectedRes < 1e-3) { isSuccess = true; diff --git a/tests/xmltester/tests/general/TestDistance.xml b/tests/xmltester/tests/general/TestDistance.xml index e6444ddac4..34e96b93e1 100644 --- a/tests/xmltester/tests/general/TestDistance.xml +++ b/tests/xmltester/tests/general/TestDistance.xml @@ -4,8 +4,8 @@ PeP - point to an empty point POINT(10 10) POINT EMPTY - 0.0 - 0.0 + NaN + NaN @@ -36,8 +36,8 @@ LL - line to empty line LINESTRING (0 0, 0 10) LINESTRING EMPTY - 0.0 - 0.0 + NaN + NaN @@ -68,8 +68,8 @@ PA - point to empty polygon POINT (240 160) POLYGON EMPTY - 0.0 - 0.0 + NaN + NaN From 6233511ce1003a7314f7eae4a24bfcff1c519e82 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Thu, 18 Dec 2025 01:53:19 +0100 Subject: [PATCH 2/2] DistanceOp::distance(): make it return inf instead of 0 if one of the geometry is empty --- capi/geos_c.h.in | 4 ++-- src/operation/distance/DistanceOp.cpp | 4 ++-- .../unit/operation/distance/DistanceOpTest.cpp | 18 +++++++++--------- tests/xmltester/XMLTester.cpp | 8 +++++++- tests/xmltester/tests/general/TestDistance.xml | 12 ++++++------ 5 files changed, 26 insertions(+), 20 deletions(-) diff --git a/capi/geos_c.h.in b/capi/geos_c.h.in index 78bb3b401f..ce40737a38 100644 --- a/capi/geos_c.h.in +++ b/capi/geos_c.h.in @@ -3631,8 +3631,8 @@ extern int GEOS_DLL GEOSGeomGetLength( * Calculate the distance between two geometries. * \param[in] g1 Input geometry * \param[in] g2 Input geometry -* \param[out] dist Pointer to be filled in with distance result. NaN is returned -* if one of the geometry is empty. +* \param[out] dist Pointer to be filled in with distance result. Positive +* infinity is returned if one of the geometry is empty. * \return 1 on success, 0 on exception. * \since 2.2 */ diff --git a/src/operation/distance/DistanceOp.cpp b/src/operation/distance/DistanceOp.cpp index cf457fbe30..2eb82eaac2 100644 --- a/src/operation/distance/DistanceOp.cpp +++ b/src/operation/distance/DistanceOp.cpp @@ -98,7 +98,7 @@ DistanceOp::DistanceOp(const Geometry& g0, const Geometry& g1, double tdist) /** * Report the distance between the closest points on the input geometries. * - * @return the distance between the geometries, or NaN if one of them is empty. + * @return the distance between the geometries, or positive infinity if one of them is empty. */ double DistanceOp::distance() @@ -112,7 +112,7 @@ DistanceOp::distance() throw IllegalArgumentException("null geometries are not supported"); } if(geom[0]->isEmpty() || geom[1]->isEmpty()) { - return std::numeric_limits::quiet_NaN(); + return std::numeric_limits::infinity(); } if(geom[0]->getGeometryTypeId() == GEOS_POINT && geom[1]->getGeometryTypeId() == GEOS_POINT) { return static_cast(geom[0])->getCoordinate()->distance(*static_cast(geom[1])->getCoordinate()); diff --git a/tests/unit/operation/distance/DistanceOpTest.cpp b/tests/unit/operation/distance/DistanceOpTest.cpp index cbd96f8060..bbb5e08504 100644 --- a/tests/unit/operation/distance/DistanceOpTest.cpp +++ b/tests/unit/operation/distance/DistanceOpTest.cpp @@ -228,7 +228,7 @@ void object::test<8> DistanceOp dist(g0.get(), g1.get()); - ensure(std::isnan(dist.distance())); + ensure(std::isinf(dist.distance())); ensure(dist.nearestPoints() == nullptr); } @@ -416,7 +416,7 @@ void object::test<16> DistanceOp dist(g0.get(), g1.get()); - ensure(std::isnan(dist.distance())); + ensure(std::isinf(dist.distance())); ensure(dist.nearestPoints() == nullptr); } @@ -497,7 +497,7 @@ void object::test<19> ensure(g1->isValid()); ensure(g2->isValid()); - ensure(std::isnan(g1->distance(g2.get()))); + ensure(std::isinf(g1->distance(g2.get()))); } // Test case reported in Shapely @@ -588,8 +588,8 @@ void object::test<23>() auto g2 = wktreader.read("LINESTRING EMPTY"); ensure(g1 != nullptr && g2 != nullptr); - ensure(std::isnan(g1->distance(g2.get()))); - ensure(std::isnan(g2->distance(g1.get()))); + ensure(std::isinf(g1->distance(g2.get()))); + ensure(std::isinf(g2->distance(g1.get()))); } template<> @@ -600,8 +600,8 @@ void object::test<24>() auto g2 = wktreader.read("LINESTRING EMPTY"); ensure(g1 != nullptr && g2 != nullptr); - ensure(std::isnan(g1->distance(g2.get()))); - ensure(std::isnan(g2->distance(g1.get()))); + ensure(std::isinf(g1->distance(g2.get()))); + ensure(std::isinf(g2->distance(g1.get()))); } // But ignore empty if there's a real distance? @@ -637,8 +637,8 @@ void object::test<27>() auto g2 = wktreader.read("GEOMETRYCOLLECTION(POINT(1 0))"); ensure(g1 != nullptr && g2 != nullptr); - ensure(std::isnan(g1->distance(g2.get()))); - ensure(std::isnan(g2->distance(g1.get()))); + ensure(std::isinf(g1->distance(g2.get()))); + ensure(std::isinf(g2->distance(g1.get()))); } diff --git a/tests/xmltester/XMLTester.cpp b/tests/xmltester/XMLTester.cpp index acdba782da..be54458132 100644 --- a/tests/xmltester/XMLTester.cpp +++ b/tests/xmltester/XMLTester.cpp @@ -896,7 +896,10 @@ Test::checkResult( double result) { char* rest = nullptr; const double expectedRes = (opResult == "NaN" || opResult == "nan") ? - std::numeric_limits::quiet_NaN() : std::strtod(opResult.c_str(), &rest); + std::numeric_limits::quiet_NaN() : + (opResult == "Inf" || opResult == "inf") ? + std::numeric_limits::infinity() : + std::strtod(opResult.c_str(), &rest); if(rest == opResult.c_str()) { throw std::runtime_error("malformed testcase: missing expected double value"); } @@ -908,6 +911,9 @@ Test::checkResult( double result) else if (std::isnan(expectedRes)) { isSuccess = std::isnan(result); } + else if( expectedRes == result ) { + isSuccess = true; + } else { if (std::abs(expectedRes - result) / expectedRes < 1e-3) { isSuccess = true; diff --git a/tests/xmltester/tests/general/TestDistance.xml b/tests/xmltester/tests/general/TestDistance.xml index 34e96b93e1..36ddb40c99 100644 --- a/tests/xmltester/tests/general/TestDistance.xml +++ b/tests/xmltester/tests/general/TestDistance.xml @@ -4,8 +4,8 @@ PeP - point to an empty point POINT(10 10) POINT EMPTY - NaN - NaN + inf + inf @@ -36,8 +36,8 @@ LL - line to empty line LINESTRING (0 0, 0 10) LINESTRING EMPTY - NaN - NaN + inf + inf @@ -68,8 +68,8 @@ PA - point to empty polygon POINT (240 160) POLYGON EMPTY - NaN - NaN + inf + inf