DistanceOp::distance(): make it return +inf instead of 0 if one of the geometry is empty (#1345)

* DistanceOp::distance(): make it return NaN instead of 0 if one of the geometry is empty

Fixes OSGeo/gdal#12978

* DistanceOp::distance(): make it return inf instead of 0 if one of the geometry is empty

Author: rouault

capi/geos_c.h.in

@@ -3405,7 +3405,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
@@ -3649,7 +3649,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. Positive
+*                  infinity is returned if one of the geometry is empty.
 * \return 1 on success, 0 on exception.
 * \since 2.2
 */
@@ -5531,7 +5532,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

src/operation/distance/DistanceOp.cpp

@@ -36,6 +36,7 @@
 #include <geos/util/IllegalArgumentException.h>
 #include <geos/util.h>
 
+#include <limits>
 #include <vector>
 #include <iostream>
 
@@ -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 positive infinity 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<double>::infinity();
     }
     if(geom[0]->getGeometryTypeId() == GEOS_POINT && geom[1]->getGeometryTypeId() == GEOS_POINT) {
         return static_cast<const Point*>(geom[0])->getCoordinate()->distance(*static_cast<const Point*>(geom[1])->getCoordinate());

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::isinf(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::isinf(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::isinf(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::isinf(g1->distance(g2.get())));
+    ensure(std::isinf(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::isinf(g1->distance(g2.get())));
+    ensure(std::isinf(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::isinf(g1->distance(g2.get())));
+    ensure(std::isinf(g2->distance(g1.get())));
 }
 
 

tests/xmltester/XMLTester.cpp

@@ -894,8 +894,12 @@ 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<double>::quiet_NaN() :
+        (opResult == "Inf" || opResult == "inf") ?
+        std::numeric_limits<double>::infinity() :
+        std::strtod(opResult.c_str(), &rest);
     if(rest == opResult.c_str()) {
         throw std::runtime_error("malformed testcase: missing expected double value");
     }
@@ -904,6 +908,12 @@ Test::checkResult( double result)
             isSuccess = true;
         }
     }
+    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;

tests/xmltester/tests/general/TestDistance.xml

@@ -4,8 +4,8 @@
   <desc>PeP - point to an empty point</desc>
   <a>    POINT(10 10)  </a>
   <b>    POINT EMPTY  </b>
-<test><op name="distance" arg1="A" arg2="B">    0.0   </op></test>
-<test><op name="distance" arg1="B" arg2="A">    0.0   </op></test>
+<test><op name="distance" arg1="A" arg2="B">    inf   </op></test>
+<test><op name="distance" arg1="B" arg2="A">    inf   </op></test>
 </case>
 
 <case>
@@ -36,8 +36,8 @@
   <desc>LL - line to empty line</desc>
   <a>    LINESTRING (0 0, 0 10)  </a>
   <b>    LINESTRING EMPTY  </b>
-<test><op name="distance" arg1="A" arg2="B">    0.0   </op></test>
-<test><op name="distance" arg1="B" arg2="A">    0.0   </op></test>
+<test><op name="distance" arg1="A" arg2="B">    inf   </op></test>
+<test><op name="distance" arg1="B" arg2="A">    inf   </op></test>
 </case>
 
 <case>
@@ -68,8 +68,8 @@
   <desc>PA - point to empty polygon</desc>
   <a>    POINT (240 160)  </a>
   <b>    POLYGON EMPTY  </b>
-<test><op name="distance" arg1="A" arg2="B" >    0.0   </op></test>
-<test><op name="distance" arg1="B" arg2="A" >    0.0   </op></test>
+<test><op name="distance" arg1="A" arg2="B" >    inf   </op></test>
+<test><op name="distance" arg1="B" arg2="A" >    inf   </op></test>
 </case>
 
 <case>