Improve scale handling for PrecisionModel

capi/geos_ts_c.cpp

@@ -2978,8 +2978,8 @@ extern "C" {
         return execute(extHandle, [&]() {
             PrecisionModel newpm;
             if(gridSize != 0) {
-                // Use negative scale to indicate you actually want a gridSize
-                double scale = -1.0 * std::abs(gridSize);
+                // Convert gridSize to scale factor
+                double scale = 1.0 / std::abs(gridSize);
                 newpm = PrecisionModel(scale);
             }
 

include/geos/geom/PrecisionModel.h

@@ -77,8 +77,6 @@ namespace geom { // geos::geom
  *
  * It is also supported to specify a precise grid size
  * by providing it as a negative scale factor.
- * This allows setting a precise grid size rather than using a fractional scale,
- * which provides more accurate and robust rounding.
  * For example, to specify rounding to the nearest 1000 use a scale factor of -1000.
  *
  * Coordinates are represented internally as Java double-precision values.
@@ -348,6 +346,11 @@ class GEOS_DLL PrecisionModel {
     void setScale(double newScale);
     // throw IllegalArgumentException
 
+    /** \brief
+     * Snaps a value to nearest integer, if within tolerance.
+     */
+    static double snapToInt(double val, double tolerance);
+
     Type modelType;
 
     /**

src/geom/PrecisionModel.cpp

@@ -27,6 +27,7 @@
 #include <string>
 #include <cmath>
 #include <iostream>
+#include <iomanip> 
 
 #ifndef GEOS_DEBUG
 #define GEOS_DEBUG 0
@@ -55,10 +56,17 @@ PrecisionModel::makePrecise(double val) const
         return static_cast<double>(floatSingleVal);
     }
     if(modelType == FIXED) {
-        if (gridSize > 0) {
+        //-- make arithmetic robust by using integral value if available
+        if (gridSize > 1) {
+//double v2 = util::round(val / gridSize) * gridSize;
+//std::cout << std::setprecision(16) << "GS[" << gridSize << "] " << val << " -> "  << v2 << std::endl;
             return util::round(val / gridSize) * gridSize;
         }
-        else {
+        //-- since grid size is <= 1, scale must be >= 1 OR 0
+        //-- if scale == 0, this is a no-op (should never happen)
+        else if (scale != 0.0) {
+//double v2 = util::round(val * scale) / scale;
+//std::cout << std::setprecision(16) << "SC[" << scale << "] " << val << " -> " << "SC " << v2 << std::endl;
             return util::round(val * scale) / scale;
         }
     }
@@ -85,7 +93,7 @@ PrecisionModel::PrecisionModel(Type nModelType)
     :
     modelType(nModelType),
     scale(1.0),
-    gridSize(0.0)
+    gridSize(1.0)
 {
 #if GEOS_DEBUG
     std::cerr << "PrecisionModel[" << this << "] ctor(Type)" << std::endl;
@@ -153,25 +161,48 @@ PrecisionModel::getMaximumSignificantDigits() const
     return maxSigDigits;
 }
 
+//-- this value is not critical, since most common usage should be VERY close to integral
+const double GRIDSIZE_INTEGER_TOLERANCE = 1e-5;
+
 /*private*/
 void
 PrecisionModel::setScale(double newScale)
 {
+    //-- should never happen, but make this a no-op in case
+    if (newScale == 0) {
+        scale = 0.0;
+        gridSize = 0.0;
+    }
     /**
     * A negative scale indicates the grid size is being set.
     * The scale is set as well, as the reciprocal.
+    * NOTE: may not need to support negative grid size now due to robust arithmetic
     */
     if (newScale < 0) {
-        gridSize = std::fabs(newScale);
-        scale = 1.0 / gridSize;
+        scale = 1.0 / std::fabs(newScale);
     }
     else {
-        scale = std::fabs(newScale);
-        /**
-        * Leave gridSize as 0, to ensure it is computed using scale
-        */
-        gridSize = 0.0;
+        scale = newScale;
+    }
+    //-- snap nearly integral scale or gridsize to exact integer
+    //-- this handles the most common case of fractional powers of ten
+    if (scale < 1) {
+        gridSize = snapToInt(1.0 / scale, GRIDSIZE_INTEGER_TOLERANCE);
     }
+    else {
+        scale = snapToInt( scale, GRIDSIZE_INTEGER_TOLERANCE);
+        gridSize = 1.0 / scale;
+    }
+}
+
+/*private*/ 
+double
+PrecisionModel::snapToInt(double val, double tolerance) {
+    double valInt = std::round(val);
+    if (std::abs(val - valInt) < tolerance) {
+        return valInt;
+    }
+    return val;
 }
 
 /*public*/

tests/unit/capi/GEOSGeom_setPrecisionTest.cpp

@@ -1,3 +1,5 @@
+//
+// Test Suite for capi::GEOSGeom_setPrecision
 
 #include "capi_test_utils.h"
 
@@ -6,8 +8,19 @@ namespace tut {
 // Test Group
 //
 
-// Common data used in test cases.
-struct test_capigeosgeomsetprecision_data : public capitest::utility {};
+// Common code used in test cases.
+struct test_capigeosgeomsetprecision_data : public capitest::utility {
+    void
+    checkPrecision(const char* wktInput, double gridSize, const char* wktExpected)
+    {
+        GEOSGeometry* input = fromWKT(wktInput);
+        GEOSGeometry* result = GEOSGeom_setPrecision(input, gridSize, 0);
+        ensure(result != nullptr);
+        ensure_geometry_equals(result, wktExpected);
+        GEOSGeom_destroy(input);
+        GEOSGeom_destroy(result);
+    }
+};
 
 typedef test_group<test_capigeosgeomsetprecision_data> group;
 typedef group::object object;
@@ -37,10 +50,9 @@ template<>
 void object::test<2>
 ()
 {
-    geom1_ = fromWKT("LINESTRING(-3 6, 9 1)");
-    geom3_ = GEOSGeom_setPrecision(geom1_, 2.0, 0);
-    ensure(geom3_ != 0);
-    ensure_geometry_equals(geom3_, "LINESTRING (-2 6, 10 2)");
+    checkPrecision("LINESTRING(-3 6, 9 1)",
+        2.0,
+        "LINESTRING (-2 6, 10 2)");
 }
 
 // See effects of precision reduction on intersection operation
@@ -123,9 +135,9 @@ template<>
 template<>
 void object::test<6> ()
 {
-    geom1_ = fromWKT("LINESTRING (0 0, 0.1 0.1)");
-    geom2_ = GEOSGeom_setPrecision(geom1_, 1.0, 0);
-    ensure_geometry_equals(geom2_, "LINESTRING EMPTY");
+    checkPrecision("LINESTRING (0 0, 0.1 0.1)",
+        1.0,
+        "LINESTRING EMPTY");
 }
 
 // Retain (or not) collapsed elements
@@ -164,41 +176,154 @@ template<>
 template<>
 void object::test<10> ()
 {
-    geom1_ = fromWKT("LINEARRING (0 0, 0.1 0, 0.1 0.1, 0 0.1, 0 0)");
-    geom2_ = GEOSGeom_setPrecision(geom1_, 1.0, 0);
-    ensure_geometry_equals(geom2_, "LINEARRING EMPTY");
+    checkPrecision("LINEARRING (0 0, 0.1 0, 0.1 0.1, 0 0.1, 0 0)",
+        1.0,
+        "LINEARRING EMPTY");
 }
 
 // Reduce polygon precision, corner case / Trac #1127
 template<>
 template<>
 void object::test<11> ()
 {
-    // POLYGON((
-    //   100 49.5, (1)
-    //   100 300,  (2)
-    //   320 60,   (3)
-    //   340 49.9, (4)
-    //   360 50.1, (5)
-    //   380 49.5, (6)
-    //   100 49.5  (7)
-    // ))
-    // * points 4 and 5 are close (less than 100.0/100) to segment (6, 7);
-    // * Y coordinates of points 4 and 5 are rounded to different values, 0 and 100 respectively;
-    // * point 4 belongs to monotone chain of size > 1 -- segments (2, 3) and (3, 4)
-    geom1_ = fromWKT("POLYGON((100 49.5, 100 300, 320 60, 340 49.9, 360 50.1, 380 49.5, 100 49.5))");
-    geom2_ = GEOSGeom_setPrecision(geom1_, 100.0, 0);
-    ensure(geom2_ != nullptr); // just check that valid geometry is constructed
+    checkPrecision("POLYGON((100 49.5, 100 300, 320 60, 340 49.9, 360 50.1, 380 49.5, 100 49.5))",
+        100.0,
+        "POLYGON ((100 300, 300 100, 300 0, 100 0, 100 300))"); 
 }
 
 template<>
 template<>
 void object::test<12>()
 {
-    geom1_ = fromWKT("POLYGON ((0 0, 0.1 0, 0.1 0.1, 0 0.1, 0 0))");
-    geom2_ = GEOSGeom_setPrecision(geom1_, 1.0, 0);
+    checkPrecision("POLYGON ((0 0, 0.1 0, 0.1 0.1, 0 0.1, 0 0))",
+        1.0,
+        "POLYGON EMPTY");
+}
+
+//-- test that a large gridsize works
+template<>
+template<>
+void object::test<13>()
+{
+    checkPrecision("LINESTRING (657035.913 6475590.114,657075.57 6475500)",
+        100,
+        "LINESTRING (657000 6475600, 657100 6475500)");
+}
+
+// Test more exact rounding for integral scale factors
+// see https://trac.osgeo.org/postgis/ticket/5520
+template<>
+template<>
+void object::test<14>()
+{
+    checkPrecision("LINESTRING (657035.913 6475590.114,657075.57 6475500)",
+        0.001,
+        "LINESTRING (657035.913 6475590.114, 657075.57 6475500)");
+}
+
+// see https://trac.osgeo.org/postgis/ticket/5425
+template<>
+template<>
+void object::test<15>()
+{
+    checkPrecision("LINESTRING(674169.89 198051.38,674197.7 198065.55,674200.36 198052.38)",
+        0.001,
+        "LINESTRING (674169.89 198051.38, 674197.7 198065.55, 674200.36 198052.38)");
+}
+
+// see https://trac.osgeo.org/postgis/ticket/3929
+template<>
+template<>
+void object::test<16>()
+{
+    checkPrecision("POINT(311.4 0)",
+        0.1,
+        "POINT(311.4 0)");
+}
+
+// see https://gis.stackexchange.com/questions/465485/postgis-reduce-precision-in-linestring
+template<>
+template<>
+void object::test<17>()
+{
+    checkPrecision("LINESTRING (16.418792399944802 54.24801559999939, 16.4176588 54.248003)",
+        0.0000001,
+        "LINESTRING (16.4187924 54.2480156, 16.4176588 54.248003)");
+}
+
+// see https://gis.stackexchange.com/questions/321814/st-snaptogrid-doesnt-work-properly-e-g-41-94186153740355-41-94186149999999
+template<>
+template<>
+void object::test<18>()
+{
+    checkPrecision("POINT (21.619820510769063 41.94186153740355)",
+        0.0000001,
+        "POINT (21.6198205 41.9418615)");
+}
+
+// see https://gis.stackexchange.com/questions/321814/st-snaptogrid-doesnt-work-properly-e-g-41-94186153740355-41-94186149999999
+template<>
+template<>
+void object::test<19>()
+{
+    checkPrecision("POINT (22.49594094391644 41.20357506925623)",
+        0.0000001, 
+        "POINT (22.4959409 41.2035751)");
+}
+
+// see https://lists.osgeo.org/pipermail/postgis-users/2006-January/010861.html
+template<>
+template<>
+void object::test<20>()
+{
+    geom1_ = fromWKT("POINT(1.23456789 9.87654321)");
+    geom2_ = GEOSGeom_setPrecision(geom1_, .000001, 0);
+    geom3_ = GEOSGeom_setPrecision(geom2_, .001, 0);
+    ensure_geometry_equals(geom3_,
+        "POINT(1.235 9.877)");
+}
+
+// see https://lists.osgeo.org/pipermail/postgis-users/2023-September/046107.html
+template<>
+template<>
+void object::test<21>()
+{
+    checkPrecision("LINESTRING(334729.13 4103548.88, 334729.12 4103548.53)",
+        0.001,
+        "LINESTRING(334729.13 4103548.88,334729.12 4103548.53)");
+}
 
-    ensure_equals(GEOSGeom_getCoordinateDimension(geom2_), 2);
+// Test multiple grid sizes
+template<>
+template<>
+void object::test<22>()
+{
+    const char* wkt = "LINESTRING(674169.89 198051.619820510769063, 674197.71234 1448065.55674200)";
+
+    checkPrecision(wkt, 0.1,       "LINESTRING (674169.9  198051.6,       674197.7     1448065.6 )");
+    checkPrecision(wkt, 0.01,      "LINESTRING (674169.89 198051.62,      674197.71    1448065.56 )");
+    checkPrecision(wkt, 0.001,     "LINESTRING (674169.89 198051.62,      674197.712   1448065.557 )");
+    checkPrecision(wkt, 0.0001,    "LINESTRING (674169.89 198051.6198,    674197.7123  1448065.5567 )");
+    checkPrecision(wkt, 0.00001,   "LINESTRING (674169.89 198051.61982,   674197.71234 1448065.55674 )");
+    checkPrecision(wkt, 0.000001,  "LINESTRING (674169.89 198051.619821,  674197.71234 1448065.556742 )");
+    checkPrecision(wkt, 0.0000001, "LINESTRING (674169.89 198051.6198205, 674197.71234 1448065.556742 )");
+
+    checkPrecision(wkt,       1, "LINESTRING ( 674170 198052,  674198 1448066)");
+    checkPrecision(wkt,      10, "LINESTRING ( 674170 198050,  674200 1448070)");
+    checkPrecision(wkt,     100, "LINESTRING ( 674200 198100,  674200 1448100)");
+    checkPrecision(wkt,    1000, "LINESTRING ( 674000 198000,  674000 1448000)");
+    checkPrecision(wkt,   10000, "LINESTRING ( 670000 200000,  670000 1450000)");
+    checkPrecision(wkt,  100000, "LINESTRING ( 700000 200000,  700000 1400000)");
+    checkPrecision(wkt, 1000000, "LINESTRING (1000000      0, 1000000 1000000)");
+}
+
+// This case with a large scale factor produced inexact rounding before code update 
+template<>
+template<>
+void object::test<23>()
+{
+    const char* wkt = "LINESTRING(674169.89 198051.619820510769063, 674197.71234 1448065.55674200)";
+    checkPrecision(wkt,  100000, "LINESTRING ( 700000 200000,  700000 1400000)");
 }
 
 } // namespace tut