diff --git a/capi/geos_ts_c.cpp b/capi/geos_ts_c.cpp index 1611118b1d..539435ffbf 100644 --- a/capi/geos_ts_c.cpp +++ b/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); } diff --git a/include/geos/geom/PrecisionModel.h b/include/geos/geom/PrecisionModel.h index e6732fa9d4..f97fa22f4d 100644 --- a/include/geos/geom/PrecisionModel.h +++ b/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; /** diff --git a/src/geom/PrecisionModel.cpp b/src/geom/PrecisionModel.cpp index 71af03a2e0..25d4085f79 100644 --- a/src/geom/PrecisionModel.cpp +++ b/src/geom/PrecisionModel.cpp @@ -27,6 +27,7 @@ #include #include #include +#include #ifndef GEOS_DEBUG #define GEOS_DEBUG 0 @@ -55,10 +56,17 @@ PrecisionModel::makePrecise(double val) const return static_cast(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*/ diff --git a/tests/unit/capi/GEOSGeom_setPrecisionTest.cpp b/tests/unit/capi/GEOSGeom_setPrecisionTest.cpp index bf792b6bff..4640b19e7e 100644 --- a/tests/unit/capi/GEOSGeom_setPrecisionTest.cpp +++ b/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 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,9 +176,9 @@ 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 @@ -174,31 +186,144 @@ 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