[geos-commits] [SCM] GEOS branch main updated. a423241d8afd797b615ed2bc8fb06c067e6aa30d
git at osgeo.org
git at osgeo.org
Wed Sep 20 15:13:58 PDT 2023
This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "GEOS".
The branch, main has been updated
via a423241d8afd797b615ed2bc8fb06c067e6aa30d (commit)
from f84966acecc6decd2aee2e7b26af83f50c6ccac3 (commit)
Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.
- Log -----------------------------------------------------------------
commit a423241d8afd797b615ed2bc8fb06c067e6aa30d
Author: Martin Davis <mtnclimb at gmail.com>
Date: Wed Sep 20 15:11:46 2023 -0700
Improve scale handling for PrecisionModel (#956)
diff --git a/capi/geos_ts_c.cpp b/capi/geos_ts_c.cpp
index 1611118b1..539435ffb 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 e6732fa9d..f97fa22f4 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 @@ private:
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 71af03a2e..25d4085f7 100644
--- a/src/geom/PrecisionModel.cpp
+++ b/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*/
diff --git a/tests/unit/capi/GEOSGeom_setPrecisionTest.cpp b/tests/unit/capi/GEOSGeom_setPrecisionTest.cpp
index bf792b6bf..4640b19e7 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<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,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");
+}
- ensure_equals(GEOSGeom_getCoordinateDimension(geom2_), 2);
+//-- 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)");
+}
+
+// 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
-----------------------------------------------------------------------
Summary of changes:
capi/geos_ts_c.cpp | 4 +-
include/geos/geom/PrecisionModel.h | 7 +-
src/geom/PrecisionModel.cpp | 51 +++++--
tests/unit/capi/GEOSGeom_setPrecisionTest.cpp | 185 +++++++++++++++++++++-----
4 files changed, 203 insertions(+), 44 deletions(-)
hooks/post-receive
--
GEOS
More information about the geos-commits
mailing list