From 10ef3b8462530ec8a6d6e21a767b7353c23b56dd Mon Sep 17 00:00:00 2001 From: Martin Davis Date: Fri, 8 Sep 2023 16:01:09 -0700 Subject: [PATCH 01/18] Improve scale handling for GEOS_setPrecision --- capi/geos_ts_c.cpp | 13 ++++++++++--- tests/unit/capi/GEOSGeom_setPrecisionTest.cpp | 13 +++++++++++++ 2 files changed, 23 insertions(+), 3 deletions(-) diff --git a/capi/geos_ts_c.cpp b/capi/geos_ts_c.cpp index 1611118b1d..bc3dd2a521 100644 --- a/capi/geos_ts_c.cpp +++ b/capi/geos_ts_c.cpp @@ -2978,9 +2978,16 @@ 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); - newpm = PrecisionModel(scale); + double scaleNom = 1.0 / gridSize; + double scaleInt = std::floor(scaleNom); + if (std::abs(scaleNom - scaleInt) < 1e-5) { + //-- if scale factor is nearly integral, use an exact integer value + newpm = PrecisionModel(scaleInt); + } + else { + //-- otherwise, just use the reciprocal of the grid size + newpm = PrecisionModel(scaleNom); + } } const PrecisionModel* pm = g->getPrecisionModel(); diff --git a/tests/unit/capi/GEOSGeom_setPrecisionTest.cpp b/tests/unit/capi/GEOSGeom_setPrecisionTest.cpp index bf792b6bff..d80c4eed79 100644 --- a/tests/unit/capi/GEOSGeom_setPrecisionTest.cpp +++ b/tests/unit/capi/GEOSGeom_setPrecisionTest.cpp @@ -201,5 +201,18 @@ void object::test<12>() ensure_equals(GEOSGeom_getCoordinateDimension(geom2_), 2); } +// Test more exact rounding for integral scale factors +// see https://trac.osgeo.org/postgis/ticket/5520 +template<> +template<> +void object::test<13>() +{ + geom1_ = fromWKT("LINESTRING (657035.913 6475590.114,657075.57 6475500)"); + geom2_ = GEOSGeom_setPrecision(geom1_, .001, 0); +//std::cout << toWKT(geom2_) << std::endl; + ensure_geometry_equals(geom2_, + "LINESTRING (657035.913 6475590.114, 657075.57 6475500)"); +} + } // namespace tut From 88de03e99b49762cc93a97b0951471dfe023abe6 Mon Sep 17 00:00:00 2001 From: Martin Davis Date: Fri, 8 Sep 2023 16:41:07 -0700 Subject: [PATCH 02/18] Refactor Grid Size integer tolerance as const --- capi/geos_ts_c.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/capi/geos_ts_c.cpp b/capi/geos_ts_c.cpp index bc3dd2a521..27ee9c0e91 100644 --- a/capi/geos_ts_c.cpp +++ b/capi/geos_ts_c.cpp @@ -2969,6 +2969,8 @@ extern "C" { }); } + const double GRIDSIZE_INTEGER_TOLERANCE = 1e-5; + Geometry* GEOSGeom_setPrecision_r(GEOSContextHandle_t extHandle, const GEOSGeometry* g, double gridSize, int flags) @@ -2980,7 +2982,7 @@ extern "C" { if(gridSize != 0) { double scaleNom = 1.0 / gridSize; double scaleInt = std::floor(scaleNom); - if (std::abs(scaleNom - scaleInt) < 1e-5) { + if (std::abs(scaleNom - scaleInt) < GRIDSIZE_INTEGER_TOLERANCE) { //-- if scale factor is nearly integral, use an exact integer value newpm = PrecisionModel(scaleInt); } From 4d2ffcbacf9efe74a7f4125a82e7bae44926f7a1 Mon Sep 17 00:00:00 2001 From: Martin Davis Date: Sat, 9 Sep 2023 08:06:47 -0700 Subject: [PATCH 03/18] Simplify logic, ensure floor works --- capi/geos_ts_c.cpp | 19 +++++++++---------- 1 file changed, 9 insertions(+), 10 deletions(-) diff --git a/capi/geos_ts_c.cpp b/capi/geos_ts_c.cpp index 27ee9c0e91..90a5637311 100644 --- a/capi/geos_ts_c.cpp +++ b/capi/geos_ts_c.cpp @@ -2979,17 +2979,16 @@ extern "C" { return execute(extHandle, [&]() { PrecisionModel newpm; - if(gridSize != 0) { - double scaleNom = 1.0 / gridSize; - double scaleInt = std::floor(scaleNom); - if (std::abs(scaleNom - scaleInt) < GRIDSIZE_INTEGER_TOLERANCE) { - //-- if scale factor is nearly integral, use an exact integer value - newpm = PrecisionModel(scaleInt); - } - else { - //-- otherwise, just use the reciprocal of the grid size - newpm = PrecisionModel(scaleNom); + if (gridSize != 0) { + //-- check for an integral scale + double scale = 1.0 / gridSize; + //-- add a small "bump" to ensure flooring works + double scaleInt = std::floor(scale + 0.1); + //-- if scale factor is essentially integral, use the exact integer value + if (std::abs(scale - scaleInt) < GRIDSIZE_INTEGER_TOLERANCE) { + scale = scaleInt; } + newpm = PrecisionModel(scale); } const PrecisionModel* pm = g->getPrecisionModel(); From 4754d776d8870e5e8779ee39dbb466853add29fb Mon Sep 17 00:00:00 2001 From: Martin Davis Date: Sat, 9 Sep 2023 17:39:53 -0700 Subject: [PATCH 04/18] Use better bump for safe floor evaluation --- capi/geos_ts_c.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/capi/geos_ts_c.cpp b/capi/geos_ts_c.cpp index 90a5637311..e0f4138a1a 100644 --- a/capi/geos_ts_c.cpp +++ b/capi/geos_ts_c.cpp @@ -2983,7 +2983,7 @@ extern "C" { //-- check for an integral scale double scale = 1.0 / gridSize; //-- add a small "bump" to ensure flooring works - double scaleInt = std::floor(scale + 0.1); + double scaleInt = std::floor(scale + 2 * GRIDSIZE_INTEGER_TOLERANCE); //-- if scale factor is essentially integral, use the exact integer value if (std::abs(scale - scaleInt) < GRIDSIZE_INTEGER_TOLERANCE) { scale = scaleInt; From 7ece4a6f1d7d48e834194862a8bf9387136d419d Mon Sep 17 00:00:00 2001 From: Martin Davis Date: Mon, 11 Sep 2023 13:32:02 -0700 Subject: [PATCH 05/18] Use std::round --- capi/geos_ts_c.cpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/capi/geos_ts_c.cpp b/capi/geos_ts_c.cpp index e0f4138a1a..04a34c28ff 100644 --- a/capi/geos_ts_c.cpp +++ b/capi/geos_ts_c.cpp @@ -2982,8 +2982,7 @@ extern "C" { if (gridSize != 0) { //-- check for an integral scale double scale = 1.0 / gridSize; - //-- add a small "bump" to ensure flooring works - double scaleInt = std::floor(scale + 2 * GRIDSIZE_INTEGER_TOLERANCE); + double scaleInt = std::round(scale); //-- if scale factor is essentially integral, use the exact integer value if (std::abs(scale - scaleInt) < GRIDSIZE_INTEGER_TOLERANCE) { scale = scaleInt; From 7173943ed356d2dbc58741dd6f31556b82693ab1 Mon Sep 17 00:00:00 2001 From: Martin Davis Date: Mon, 11 Sep 2023 15:01:28 -0700 Subject: [PATCH 06/18] Only check integral scale facter for gridsize lt 1 --- capi/geos_ts_c.cpp | 12 +++++++----- tests/unit/capi/GEOSGeom_setPrecisionTest.cpp | 12 ++++++++++++ 2 files changed, 19 insertions(+), 5 deletions(-) diff --git a/capi/geos_ts_c.cpp b/capi/geos_ts_c.cpp index 04a34c28ff..da11e3690c 100644 --- a/capi/geos_ts_c.cpp +++ b/capi/geos_ts_c.cpp @@ -2980,12 +2980,14 @@ extern "C" { return execute(extHandle, [&]() { PrecisionModel newpm; if (gridSize != 0) { - //-- check for an integral scale double scale = 1.0 / gridSize; - double scaleInt = std::round(scale); - //-- if scale factor is essentially integral, use the exact integer value - if (std::abs(scale - scaleInt) < GRIDSIZE_INTEGER_TOLERANCE) { - scale = scaleInt; + //-- check if the gridsize corresponds to an integral scale + if (gridSize < 1) { + double scaleInt = std::round(scale); + //-- if scale factor is essentially integral, use the exact integer value + if (std::abs(scale - scaleInt) < GRIDSIZE_INTEGER_TOLERANCE) { + scale = scaleInt; + } } newpm = PrecisionModel(scale); } diff --git a/tests/unit/capi/GEOSGeom_setPrecisionTest.cpp b/tests/unit/capi/GEOSGeom_setPrecisionTest.cpp index d80c4eed79..cf0d00cbee 100644 --- a/tests/unit/capi/GEOSGeom_setPrecisionTest.cpp +++ b/tests/unit/capi/GEOSGeom_setPrecisionTest.cpp @@ -214,5 +214,17 @@ void object::test<13>() "LINESTRING (657035.913 6475590.114, 657075.57 6475500)"); } +//-- test that a large gridsize works +template<> +template<> +void object::test<14>() +{ + geom1_ = fromWKT("LINESTRING (657035.913 6475590.114,657075.57 6475500)"); + geom2_ = GEOSGeom_setPrecision(geom1_, 100, 0); +//std::cout << toWKT(geom2_) << std::endl; + ensure_geometry_equals(geom2_, + "LINESTRING (657000 6475600, 657100 6475500)"); +} + } // namespace tut From d43f7e09f070c73a682fda9e346c4944c61cfb93 Mon Sep 17 00:00:00 2001 From: Martin Davis Date: Tue, 12 Sep 2023 10:01:15 -0700 Subject: [PATCH 07/18] Add test from PostGIS SnapToGrid issue --- tests/unit/capi/GEOSGeom_setPrecisionTest.cpp | 26 ++++++++++++++----- 1 file changed, 20 insertions(+), 6 deletions(-) diff --git a/tests/unit/capi/GEOSGeom_setPrecisionTest.cpp b/tests/unit/capi/GEOSGeom_setPrecisionTest.cpp index cf0d00cbee..efa4b06c50 100644 --- a/tests/unit/capi/GEOSGeom_setPrecisionTest.cpp +++ b/tests/unit/capi/GEOSGeom_setPrecisionTest.cpp @@ -201,11 +201,23 @@ void object::test<12>() ensure_equals(GEOSGeom_getCoordinateDimension(geom2_), 2); } +//-- test that a large gridsize works +template<> +template<> +void object::test<13>() +{ + geom1_ = fromWKT("LINESTRING (657035.913 6475590.114,657075.57 6475500)"); + geom2_ = GEOSGeom_setPrecision(geom1_, 100, 0); +//std::cout << toWKT(geom2_) << std::endl; + ensure_geometry_equals(geom2_, + "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<13>() +void object::test<14>() { geom1_ = fromWKT("LINESTRING (657035.913 6475590.114,657075.57 6475500)"); geom2_ = GEOSGeom_setPrecision(geom1_, .001, 0); @@ -214,17 +226,19 @@ void object::test<13>() "LINESTRING (657035.913 6475590.114, 657075.57 6475500)"); } -//-- test that a large gridsize works +// see https://trac.osgeo.org/postgis/ticket/5425 template<> template<> -void object::test<14>() +void object::test<15>() { - geom1_ = fromWKT("LINESTRING (657035.913 6475590.114,657075.57 6475500)"); - geom2_ = GEOSGeom_setPrecision(geom1_, 100, 0); + geom1_ = fromWKT("LINESTRING(674169.89 198051.38,674197.7 198065.55,674200.36 198052.38)"); + geom2_ = GEOSGeom_setPrecision(geom1_, .001, 0); //std::cout << toWKT(geom2_) << std::endl; ensure_geometry_equals(geom2_, - "LINESTRING (657000 6475600, 657100 6475500)"); + "LINESTRING (674169.89 198051.38, 674197.7 198065.55, 674200.36 198052.38)"); } + + } // namespace tut From 6c7c52c6bc74616df24d93749aa1b4b83e6a4243 Mon Sep 17 00:00:00 2001 From: Martin Davis Date: Tue, 12 Sep 2023 11:40:14 -0700 Subject: [PATCH 08/18] Add SnapToGrid test case --- tests/unit/capi/GEOSGeom_setPrecisionTest.cpp | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/tests/unit/capi/GEOSGeom_setPrecisionTest.cpp b/tests/unit/capi/GEOSGeom_setPrecisionTest.cpp index efa4b06c50..8588642dc8 100644 --- a/tests/unit/capi/GEOSGeom_setPrecisionTest.cpp +++ b/tests/unit/capi/GEOSGeom_setPrecisionTest.cpp @@ -238,7 +238,17 @@ void object::test<15>() "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>() +{ + geom1_ = fromWKT("POINT(311.4 0)"); + geom2_ = GEOSGeom_setPrecision(geom1_, .1, 0); +//std::cout << toWKT(geom2_) << std::endl; + ensure_geometry_equals(geom2_, + "POINT(311.4 0)"); +} } // namespace tut From 0662640c02b114fc35e33358a7f5c3aaeefd658a Mon Sep 17 00:00:00 2001 From: Martin Davis Date: Tue, 12 Sep 2023 15:39:25 -0700 Subject: [PATCH 09/18] Add more tests --- tests/unit/capi/GEOSGeom_setPrecisionTest.cpp | 49 +++++++++++++++++-- 1 file changed, 45 insertions(+), 4 deletions(-) diff --git a/tests/unit/capi/GEOSGeom_setPrecisionTest.cpp b/tests/unit/capi/GEOSGeom_setPrecisionTest.cpp index 8588642dc8..f700c274de 100644 --- a/tests/unit/capi/GEOSGeom_setPrecisionTest.cpp +++ b/tests/unit/capi/GEOSGeom_setPrecisionTest.cpp @@ -208,7 +208,6 @@ void object::test<13>() { geom1_ = fromWKT("LINESTRING (657035.913 6475590.114,657075.57 6475500)"); geom2_ = GEOSGeom_setPrecision(geom1_, 100, 0); -//std::cout << toWKT(geom2_) << std::endl; ensure_geometry_equals(geom2_, "LINESTRING (657000 6475600, 657100 6475500)"); } @@ -221,7 +220,6 @@ void object::test<14>() { geom1_ = fromWKT("LINESTRING (657035.913 6475590.114,657075.57 6475500)"); geom2_ = GEOSGeom_setPrecision(geom1_, .001, 0); -//std::cout << toWKT(geom2_) << std::endl; ensure_geometry_equals(geom2_, "LINESTRING (657035.913 6475590.114, 657075.57 6475500)"); } @@ -233,7 +231,6 @@ void object::test<15>() { geom1_ = fromWKT("LINESTRING(674169.89 198051.38,674197.7 198065.55,674200.36 198052.38)"); geom2_ = GEOSGeom_setPrecision(geom1_, .001, 0); -//std::cout << toWKT(geom2_) << std::endl; ensure_geometry_equals(geom2_, "LINESTRING (674169.89 198051.38, 674197.7 198065.55, 674200.36 198052.38)"); } @@ -245,10 +242,54 @@ void object::test<16>() { geom1_ = fromWKT("POINT(311.4 0)"); geom2_ = GEOSGeom_setPrecision(geom1_, .1, 0); -//std::cout << toWKT(geom2_) << std::endl; ensure_geometry_equals(geom2_, "POINT(311.4 0)"); } +// see https://gis.stackexchange.com/questions/465485/postgis-reduce-precision-in-linestring +template<> +template<> +void object::test<17>() +{ + geom1_ = fromWKT("LINESTRING (16.418792399944802 54.24801559999939, 16.4176588 54.248003)"); + geom2_ = GEOSGeom_setPrecision(geom1_, .0000001, 0); + ensure_geometry_equals(geom2_, + "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>() +{ + geom1_ = fromWKT("POINT (21.619820510769063 41.94186153740355)"); + geom2_ = GEOSGeom_setPrecision(geom1_, .0000001, 0); + ensure_geometry_equals(geom2_, + "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>() +{ + geom1_ = fromWKT("POINT (22.49594094391644 41.20357506925623)"); + geom2_ = GEOSGeom_setPrecision(geom1_, .0000001, 0); + ensure_geometry_equals(geom2_, + "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)"); +} + } // namespace tut From f97c617cb5b4ab0578b5d0856ab0539381c59d7a Mon Sep 17 00:00:00 2001 From: Martin Davis Date: Tue, 12 Sep 2023 15:43:33 -0700 Subject: [PATCH 10/18] Add test case --- tests/unit/capi/GEOSGeom_setPrecisionTest.cpp | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/tests/unit/capi/GEOSGeom_setPrecisionTest.cpp b/tests/unit/capi/GEOSGeom_setPrecisionTest.cpp index f700c274de..9e908fbec7 100644 --- a/tests/unit/capi/GEOSGeom_setPrecisionTest.cpp +++ b/tests/unit/capi/GEOSGeom_setPrecisionTest.cpp @@ -291,5 +291,16 @@ void object::test<20>() "POINT(1.235 9.877)"); } +// see https://lists.osgeo.org/pipermail/postgis-users/2023-September/046107.html +template<> +template<> +void object::test<21>() +{ + geom1_ = fromWKT("LINESTRING(334729.13 4103548.88, 334729.12 4103548.53)"); + geom2_ = GEOSGeom_setPrecision(geom1_, .001, 0); + ensure_geometry_equals(geom2_, + "LINESTRING(334729.13 4103548.88,334729.12 4103548.53)"); +} + } // namespace tut From a4f4c612a937317236ccd69dae4c40de17c2b0bf Mon Sep 17 00:00:00 2001 From: Martin Davis Date: Tue, 12 Sep 2023 18:38:54 -0700 Subject: [PATCH 11/18] Add utility test function --- tests/unit/capi/GEOSGeom_setPrecisionTest.cpp | 99 +++++++++---------- 1 file changed, 44 insertions(+), 55 deletions(-) diff --git a/tests/unit/capi/GEOSGeom_setPrecisionTest.cpp b/tests/unit/capi/GEOSGeom_setPrecisionTest.cpp index 9e908fbec7..37d66f0831 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,17 @@ 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) + { + geom1_ = fromWKT(wktInput); + geom2_ = GEOSGeom_setPrecision(geom1_, gridSize, 0); + ensure(geom2_ != nullptr); + ensure_geometry_equals(geom2_, wktExpected); + } +}; typedef test_group group; typedef group::object object; @@ -37,10 +48,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 +133,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 +174,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 +184,18 @@ 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); - - ensure_equals(GEOSGeom_getCoordinateDimension(geom2_), 2); + 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 @@ -206,9 +203,8 @@ template<> template<> void object::test<13>() { - geom1_ = fromWKT("LINESTRING (657035.913 6475590.114,657075.57 6475500)"); - geom2_ = GEOSGeom_setPrecision(geom1_, 100, 0); - ensure_geometry_equals(geom2_, + checkPrecision("LINESTRING (657035.913 6475590.114,657075.57 6475500)", + 100, "LINESTRING (657000 6475600, 657100 6475500)"); } @@ -218,9 +214,8 @@ template<> template<> void object::test<14>() { - geom1_ = fromWKT("LINESTRING (657035.913 6475590.114,657075.57 6475500)"); - geom2_ = GEOSGeom_setPrecision(geom1_, .001, 0); - ensure_geometry_equals(geom2_, + checkPrecision("LINESTRING (657035.913 6475590.114,657075.57 6475500)", + 0.001, "LINESTRING (657035.913 6475590.114, 657075.57 6475500)"); } @@ -229,9 +224,8 @@ template<> template<> void object::test<15>() { - geom1_ = fromWKT("LINESTRING(674169.89 198051.38,674197.7 198065.55,674200.36 198052.38)"); - geom2_ = GEOSGeom_setPrecision(geom1_, .001, 0); - ensure_geometry_equals(geom2_, + 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)"); } @@ -240,9 +234,8 @@ template<> template<> void object::test<16>() { - geom1_ = fromWKT("POINT(311.4 0)"); - geom2_ = GEOSGeom_setPrecision(geom1_, .1, 0); - ensure_geometry_equals(geom2_, + checkPrecision("POINT(311.4 0)", + 0.1, "POINT(311.4 0)"); } @@ -251,9 +244,8 @@ template<> template<> void object::test<17>() { - geom1_ = fromWKT("LINESTRING (16.418792399944802 54.24801559999939, 16.4176588 54.248003)"); - geom2_ = GEOSGeom_setPrecision(geom1_, .0000001, 0); - ensure_geometry_equals(geom2_, + checkPrecision("LINESTRING (16.418792399944802 54.24801559999939, 16.4176588 54.248003)", + 0.0000001, "LINESTRING (16.4187924 54.2480156, 16.4176588 54.248003)"); } @@ -262,9 +254,8 @@ template<> template<> void object::test<18>() { - geom1_ = fromWKT("POINT (21.619820510769063 41.94186153740355)"); - geom2_ = GEOSGeom_setPrecision(geom1_, .0000001, 0); - ensure_geometry_equals(geom2_, + checkPrecision("POINT (21.619820510769063 41.94186153740355)", + 0.0000001, "POINT (21.6198205 41.9418615)"); } @@ -273,9 +264,8 @@ template<> template<> void object::test<19>() { - geom1_ = fromWKT("POINT (22.49594094391644 41.20357506925623)"); - geom2_ = GEOSGeom_setPrecision(geom1_, .0000001, 0); - ensure_geometry_equals(geom2_, + checkPrecision("POINT (22.49594094391644 41.20357506925623)", + 0.0000001, "POINT (22.4959409 41.2035751)"); } @@ -296,9 +286,8 @@ template<> template<> void object::test<21>() { - geom1_ = fromWKT("LINESTRING(334729.13 4103548.88, 334729.12 4103548.53)"); - geom2_ = GEOSGeom_setPrecision(geom1_, .001, 0); - ensure_geometry_equals(geom2_, + checkPrecision("LINESTRING(334729.13 4103548.88, 334729.12 4103548.53)", + 0.001, "LINESTRING(334729.13 4103548.88,334729.12 4103548.53)"); } From eaf387ae15cdd6f3ecc471c775f2315fa78735fa Mon Sep 17 00:00:00 2001 From: Martin Davis Date: Tue, 12 Sep 2023 20:24:16 -0700 Subject: [PATCH 12/18] Add test with multiple scales --- tests/unit/capi/GEOSGeom_setPrecisionTest.cpp | 21 +++++++++++++++++++ 1 file changed, 21 insertions(+) diff --git a/tests/unit/capi/GEOSGeom_setPrecisionTest.cpp b/tests/unit/capi/GEOSGeom_setPrecisionTest.cpp index 37d66f0831..e75aa7862e 100644 --- a/tests/unit/capi/GEOSGeom_setPrecisionTest.cpp +++ b/tests/unit/capi/GEOSGeom_setPrecisionTest.cpp @@ -291,5 +291,26 @@ void object::test<21>() "LINESTRING(334729.13 4103548.88,334729.12 4103548.53)"); } +// 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, 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)"); + // this fails for some reason although WKT looks identical + //checkPrecision(wkt, 100000, "LINESTRING (700000 200000, 700000 1400000)"); + checkPrecision(wkt, 1000000, "LINESTRING (1000000 0, 1000000 1000000)"); + +} } // namespace tut From 7cdbd3581b46943b9849a91476b7ce2003ab144a Mon Sep 17 00:00:00 2001 From: Martin Davis Date: Wed, 13 Sep 2023 18:29:14 -0700 Subject: [PATCH 13/18] Move rounding logic into PrecisionModel --- capi/geos_ts_c.cpp | 15 ++----- include/geos/geom/PrecisionModel.h | 5 +++ src/geom/PrecisionModel.cpp | 41 +++++++++++++++---- tests/unit/capi/GEOSGeom_setPrecisionTest.cpp | 16 +++++--- 4 files changed, 51 insertions(+), 26 deletions(-) diff --git a/capi/geos_ts_c.cpp b/capi/geos_ts_c.cpp index da11e3690c..1f02e7adaf 100644 --- a/capi/geos_ts_c.cpp +++ b/capi/geos_ts_c.cpp @@ -2969,7 +2969,7 @@ extern "C" { }); } - const double GRIDSIZE_INTEGER_TOLERANCE = 1e-5; + //const double GRIDSIZE_INTEGER_TOLERANCE = 1e-5; Geometry* GEOSGeom_setPrecision_r(GEOSContextHandle_t extHandle, const GEOSGeometry* g, @@ -2979,16 +2979,9 @@ extern "C" { return execute(extHandle, [&]() { PrecisionModel newpm; - if (gridSize != 0) { - double scale = 1.0 / gridSize; - //-- check if the gridsize corresponds to an integral scale - if (gridSize < 1) { - double scaleInt = std::round(scale); - //-- if scale factor is essentially integral, use the exact integer value - if (std::abs(scale - scaleInt) < GRIDSIZE_INTEGER_TOLERANCE) { - scale = scaleInt; - } - } + if(gridSize != 0) { + // Use negative scale to indicate argument is gridSize + 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..50151e199d 100644 --- a/include/geos/geom/PrecisionModel.h +++ b/include/geos/geom/PrecisionModel.h @@ -348,6 +348,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..815f2b0e2d 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,15 @@ 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 { +//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 +91,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,6 +159,9 @@ 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) @@ -162,18 +171,32 @@ PrecisionModel::setScale(double newScale) * The scale is set as well, as the reciprocal. */ if (newScale < 0) { - gridSize = std::fabs(newScale); - scale = 1.0 / gridSize; + scale = 1.0 / std::fabs(newScale); + } + else { + 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 = std::fabs(newScale); - /** - * Leave gridSize as 0, to ensure it is computed using scale - */ - gridSize = 0.0; + 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*/ double PrecisionModel::getOffsetX() const diff --git a/tests/unit/capi/GEOSGeom_setPrecisionTest.cpp b/tests/unit/capi/GEOSGeom_setPrecisionTest.cpp index e75aa7862e..43c7ea7b5a 100644 --- a/tests/unit/capi/GEOSGeom_setPrecisionTest.cpp +++ b/tests/unit/capi/GEOSGeom_setPrecisionTest.cpp @@ -298,19 +298,23 @@ 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.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)"); - // this fails for some reason although WKT looks identical - //checkPrecision(wkt, 100000, "LINESTRING (700000 200000, 700000 1400000)"); + // this fails due to rounding error if computed directly with scale factor + checkPrecision(wkt, 100000, "LINESTRING ( 700000 200000, 700000 1400000)"); checkPrecision(wkt, 1000000, "LINESTRING (1000000 0, 1000000 1000000)"); - } + } // namespace tut From 1b3f1af65b1465c1dc8f18e7575a712becf71343 Mon Sep 17 00:00:00 2001 From: Martin Davis Date: Wed, 13 Sep 2023 18:49:43 -0700 Subject: [PATCH 14/18] Avoid use of negative grid size for PM --- capi/geos_ts_c.cpp | 4 ++-- src/geom/PrecisionModel.cpp | 1 + 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/capi/geos_ts_c.cpp b/capi/geos_ts_c.cpp index 1f02e7adaf..510c7d8deb 100644 --- a/capi/geos_ts_c.cpp +++ b/capi/geos_ts_c.cpp @@ -2980,8 +2980,8 @@ extern "C" { return execute(extHandle, [&]() { PrecisionModel newpm; if(gridSize != 0) { - // Use negative scale to indicate argument is 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/src/geom/PrecisionModel.cpp b/src/geom/PrecisionModel.cpp index 815f2b0e2d..60743c86d6 100644 --- a/src/geom/PrecisionModel.cpp +++ b/src/geom/PrecisionModel.cpp @@ -169,6 +169,7 @@ PrecisionModel::setScale(double newScale) /** * 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) { scale = 1.0 / std::fabs(newScale); From 7d0214d7799b445d8548c7e973fc80ee9c0f590e Mon Sep 17 00:00:00 2001 From: Martin Davis Date: Wed, 13 Sep 2023 22:21:21 -0700 Subject: [PATCH 15/18] Make helper function memsafe --- tests/unit/capi/GEOSGeom_setPrecisionTest.cpp | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/tests/unit/capi/GEOSGeom_setPrecisionTest.cpp b/tests/unit/capi/GEOSGeom_setPrecisionTest.cpp index 43c7ea7b5a..4f58117fa1 100644 --- a/tests/unit/capi/GEOSGeom_setPrecisionTest.cpp +++ b/tests/unit/capi/GEOSGeom_setPrecisionTest.cpp @@ -13,10 +13,12 @@ struct test_capigeosgeomsetprecision_data : public capitest::utility { void checkPrecision(const char* wktInput, double gridSize, const char* wktExpected) { - geom1_ = fromWKT(wktInput); - geom2_ = GEOSGeom_setPrecision(geom1_, gridSize, 0); - ensure(geom2_ != nullptr); - ensure_geometry_equals(geom2_, 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); } }; From 7d5a5f883a9d25073d64b6082a90207e6bc21002 Mon Sep 17 00:00:00 2001 From: Martin Davis Date: Thu, 14 Sep 2023 10:34:55 -0700 Subject: [PATCH 16/18] Remove dead code --- capi/geos_ts_c.cpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/capi/geos_ts_c.cpp b/capi/geos_ts_c.cpp index 510c7d8deb..539435ffbf 100644 --- a/capi/geos_ts_c.cpp +++ b/capi/geos_ts_c.cpp @@ -2969,8 +2969,6 @@ extern "C" { }); } - //const double GRIDSIZE_INTEGER_TOLERANCE = 1e-5; - Geometry* GEOSGeom_setPrecision_r(GEOSContextHandle_t extHandle, const GEOSGeometry* g, double gridSize, int flags) From 778c028b44d9891133104a8f375d3f2d819cb1dd Mon Sep 17 00:00:00 2001 From: Martin Davis Date: Thu, 14 Sep 2023 10:46:59 -0700 Subject: [PATCH 17/18] Handle scale = 0 situation --- include/geos/geom/PrecisionModel.h | 2 -- src/geom/PrecisionModel.cpp | 9 ++++++++- 2 files changed, 8 insertions(+), 3 deletions(-) diff --git a/include/geos/geom/PrecisionModel.h b/include/geos/geom/PrecisionModel.h index 50151e199d..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. diff --git a/src/geom/PrecisionModel.cpp b/src/geom/PrecisionModel.cpp index 60743c86d6..25d4085f79 100644 --- a/src/geom/PrecisionModel.cpp +++ b/src/geom/PrecisionModel.cpp @@ -62,7 +62,9 @@ PrecisionModel::makePrecise(double val) const //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; @@ -166,6 +168,11 @@ const double GRIDSIZE_INTEGER_TOLERANCE = 1e-5; 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. From 1bede1e05e7e256ffcd22fd04f4f4fff1596ec06 Mon Sep 17 00:00:00 2001 From: Martin Davis Date: Thu, 14 Sep 2023 10:58:58 -0700 Subject: [PATCH 18/18] Add test case for large scale factor --- tests/unit/capi/GEOSGeom_setPrecisionTest.cpp | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/tests/unit/capi/GEOSGeom_setPrecisionTest.cpp b/tests/unit/capi/GEOSGeom_setPrecisionTest.cpp index 4f58117fa1..4640b19e7e 100644 --- a/tests/unit/capi/GEOSGeom_setPrecisionTest.cpp +++ b/tests/unit/capi/GEOSGeom_setPrecisionTest.cpp @@ -313,10 +313,18 @@ void object::test<22>() checkPrecision(wkt, 100, "LINESTRING ( 674200 198100, 674200 1448100)"); checkPrecision(wkt, 1000, "LINESTRING ( 674000 198000, 674000 1448000)"); checkPrecision(wkt, 10000, "LINESTRING ( 670000 200000, 670000 1450000)"); - // this fails due to rounding error if computed directly with scale factor 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