From becbb5c0561f85ab9fe17c5caf2ecad24fa830cb Mon Sep 17 00:00:00 2001 From: Tom Russell Date: Fri, 24 May 2024 14:26:01 +0100 Subject: [PATCH 1/4] Test for pN and pE length consistently when splitting line The specific case added as a test was triggering the condition where pE != pN but pE.length() == pN.length() and falling through to the "Unexpected points..." error. If we consistently check lengths through all of these cases, this is fixed. --- extension/src/grid.hpp | 2 +- tests/core/test_core.py | 30 ++++++++++++++++++++++++++++++ 2 files changed, 31 insertions(+), 1 deletion(-) diff --git a/extension/src/grid.hpp b/extension/src/grid.hpp index 77a4729..74168ab 100644 --- a/extension/src/grid.hpp +++ b/extension/src/grid.hpp @@ -155,7 +155,7 @@ struct Grid { // graticule crossings. If both crossing points overlap then we // update both and it doesn't matter which one we add to the // vector of crossings. - if (pE == pN) { + if (pE.length() == pN.length()) { crossings.push_back(line.start + pN); // Update the distance to the next graticule. dE += step_x; diff --git a/tests/core/test_core.py b/tests/core/test_core.py index c50b916..5789871 100644 --- a/tests/core/test_core.py +++ b/tests/core/test_core.py @@ -1,4 +1,5 @@ import pytest +import shapely.wkt import snail.core.intersections from shapely.geometry import LineString, Polygon @@ -38,6 +39,35 @@ def test_linestring_splitting(test_linestring, expected): assert split.equals_exact(expected_split, 1e-7) +def test_linestring_splitting_issue_61(): + ncols = 120163 + nrows = 259542 + transform = (5.0, 0.0, 54675.0, 0.0, -5.0, 1217320.0) + + # reduced test case to single straight-line segment + test_linestring = LineString( + [(415805.57, 430046.95), (415800.0, 430015.0)] + ) + expected = [ + LineString([(415805.57, 430046.95), (415805.23004, 430045.0)]), + LineString([(415805.23004, 430045.0), (415805, 430043.68043)]), + LineString([(415805.0, 430043.68043), (415804.35837, 430040.0)]), + LineString([(415804.35837, 430040.0), (415803.48669, 430035.0)]), + LineString([(415803.48669, 430035.0), (415802.61502, 430030.0)]), + LineString([(415802.61502, 430030.0), (415801.74334, 430025.0)]), + LineString([(415801.74334, 430025.0), (415800.87167, 430020.0)]), + LineString([(415800.87167, 430020.0), (415800.0, 430015.0)]), + ] + actual = snail.core.intersections.split_linestring( + test_linestring, nrows, ncols, transform + ) + for split, expected_split in zip(actual, expected): + if not split.equals_exact(expected_split, 1e-5): + assert ( + False + ), f"Expected split coordinates to match, got {split}, {expected_split}" + + @pytest.mark.parametrize( "test_linestring,expected", [ From a5aebfeadb74e8c00c1acd1c2344f70d8f416fea Mon Sep 17 00:00:00 2001 From: Tom Russell Date: Fri, 24 May 2024 17:52:27 +0100 Subject: [PATCH 2/4] Use std::hypot for 2D length calculation Should be more robust to numerical edge cases (0/inf/nan) and under/overflow. --- extension/src/geometry.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/extension/src/geometry.hpp b/extension/src/geometry.hpp index 3f4fd2e..443549b 100644 --- a/extension/src/geometry.hpp +++ b/extension/src/geometry.hpp @@ -33,7 +33,7 @@ struct Coord { : x(std::get<0>(xy)), y(std::get<1>(xy)) {} // Helper function to calculate the length of a Coord - inline double length(void) const { return sqrt(x * x + y * y); } + inline double length(void) const { return std::hypot(x, y); } }; /// A templated 2D line representation @@ -51,7 +51,7 @@ struct Line { inline double length(void) const { double dx = end.x - start.x; double dy = end.y - start.y; - return sqrt(dx * dx + dy * dy); + return std::hypot(dx, dy); } /// Calculate the bearing of a line. inline double bearing(void) const { From d30e186aa03dc6c2df23f93f4c38a9d1ea947694 Mon Sep 17 00:00:00 2001 From: Tom Russell Date: Fri, 24 May 2024 17:52:54 +0100 Subject: [PATCH 3/4] Drop commented debug lines --- extension/src/grid.hpp | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/extension/src/grid.hpp b/extension/src/grid.hpp index 74168ab..3532632 100644 --- a/extension/src/grid.hpp +++ b/extension/src/grid.hpp @@ -29,14 +29,14 @@ struct Grid { std::vector data; Grid() - : ncols{0}, nrows{0}, - grid_to_world{Affine()}, data{std::vector()} { + : ncols{0}, nrows{0}, grid_to_world{Affine()}, + data{std::vector()} { world_to_grid = ~grid_to_world; }; Grid(size_t ncols, size_t nrows, Affine grid_to_world) - : ncols{ncols}, nrows{nrows}, - grid_to_world{grid_to_world}, data{std::vector()} { + : ncols{ncols}, nrows{nrows}, grid_to_world{grid_to_world}, + data{std::vector()} { world_to_grid = ~grid_to_world; }; @@ -125,8 +125,6 @@ struct Grid { // Consider each possible crossing point along horizontal line if (rise == 0) { while (pE.length() <= length) { - // std::cout << " pE(" << pE.x << "," << pE.y << ") length " << - // pE.length() << "\n"; // Add the closest crossing point crossings.push_back(line.start + pE); From fd501cb60577d1bfe65faf6081dfbdcdd9e5b26b Mon Sep 17 00:00:00 2001 From: Tom Russell Date: Fri, 24 May 2024 17:56:32 +0100 Subject: [PATCH 4/4] Use cell size as a reference value in "almost-equals" comparison Issue 53 was failing with an unreasonably precise comparison of values (point, gridline) near zero. Fixed by including the grid cell size as part of the calculation of a reasonably scaled epsilon, in general the utility almost_equals now also takes an argument for an indicative order-of-magnitude reference value. --- extension/src/operations.cpp | 19 ++++++++++++++----- extension/src/utils.hpp | 19 ++++++++++++++----- tests/core/test_core.py | 24 +++++++++++++++++++----- 3 files changed, 47 insertions(+), 15 deletions(-) diff --git a/extension/src/operations.cpp b/extension/src/operations.cpp index be8681b..2415269 100644 --- a/extension/src/operations.cpp +++ b/extension/src/operations.cpp @@ -61,13 +61,14 @@ findIntersectionsLineString(geometry::LineString linestring, return (allsplits); } -bool isOnGridLine(geometry::Coord point, Direction direction, double level) { +bool isOnGridLine(geometry::Coord point, Direction direction, double level, + double cellSize) { switch (direction) { case Direction::horizontal: - return snail::utils::almost_equal(point.y, level, 2); + return snail::utils::almost_equal(point.y, level, cellSize); return (point.y == level); case Direction::vertical: - return snail::utils::almost_equal(point.x, level, 2); + return snail::utils::almost_equal(point.x, level, cellSize); return (point.x == level); default: return false; @@ -123,6 +124,8 @@ std::vector splitAlongGridlines(linestr exterior_crossings, std::vector crossings_on_gridline; std::vector gridline_splits; + double cell_size = std::max(grid.grid_to_world.a, grid.grid_to_world.e); + for (int level = min_level; level <= max_level; level++) { // find level value in coordinates double level_value = gridCoordinate(level, direction, grid); @@ -147,7 +150,7 @@ std::vector splitAlongGridlines(linestr exterior_crossings, : (curr + 1); // include if on the current line and prev/next are on opposite sides - if (isOnGridLine(*curr, direction, level_value) && + if (isOnGridLine(*curr, direction, level_value, cell_size) && crossesGridLine(*prev, *next, direction, level_value)) { crossings_on_gridline.push_back(*curr); } @@ -167,7 +170,13 @@ std::vector splitAlongGridlines(linestr exterior_crossings, }); if (crossings_on_gridline.size() % 2 != 0) { - utils::Exception("Expected even number of crossings on gridline."); + std::ostringstream msg; + msg << "Expected even number of crossings on gridline.\n"; + msg << " Found crossings on gridline:\n"; + for (auto c : crossings_on_gridline) { + msg << " c(" << c.x << "," << c.y << ")\n"; + } + utils::Exception(msg.str()); break; } diff --git a/extension/src/utils.hpp b/extension/src/utils.hpp index e16afb0..b2ad23d 100644 --- a/extension/src/utils.hpp +++ b/extension/src/utils.hpp @@ -38,13 +38,22 @@ class Exception : std::exception { /// https://en.cppreference.com/w/cpp/types/numeric_limits/epsilon template typename std::enable_if::is_integer, bool>::type -almost_equal(T x, T y, int ulp) { +almost_equal(T x, T y, T reference_value) { // the machine epsilon has to be scaled to the magnitude of the values used // and multiplied by the desired precision in ULPs (units in the last place) - return std::fabs(x - y) <= - std::numeric_limits::epsilon() * std::fabs(x + y) * ulp - // unless the result is subnormal - || std::fabs(x - y) < std::numeric_limits::min(); + // unless the result is subnormal + int ulp = 3; // magic three value for small number of units in the last place + auto abs_diff = std::fabs(x - y); + auto epsilon = std::numeric_limits::epsilon(); + // add a reference value for indicative scale, in case our x and y are + // unreasonably near zero + auto abs_total = (std::fabs(x) + std::fabs(y) + std::fabs(reference_value)); + auto scaled_epsilon = epsilon * abs_total * T(ulp); + + bool check_normal = abs_diff <= scaled_epsilon; + bool check_subnormal = abs_diff < std::numeric_limits::min(); + + return check_normal || check_subnormal; } } // namespace utils diff --git a/tests/core/test_core.py b/tests/core/test_core.py index 5789871..a550301 100644 --- a/tests/core/test_core.py +++ b/tests/core/test_core.py @@ -88,19 +88,33 @@ def test_get_cell_indices(test_linestring, expected): assert cell_indices == expected -@pytest.mark.xfail -def test_split_polygons(): +def test_split_polygons_issue_53(): + # reduced test case + polygon = Polygon( + ( + [-0.1, 2], + [-0.1, 1], + [0.1, 1], + ) + ) + + snail.core.intersections.split_polygon( + polygon, + 2, + 2, + (10.0, 0.0, -100.0, 0.0, -10.0, 100.0), + ) + bad_poly = Polygon( ( [-0.0062485600499826, 51.61041647955], [-0.0062485600499826, 51.602083146149994], [0.0020847733500204, 51.602083146149994], - # [0.0020847733500204, 51.61041647955], - # [-0.0062485600499826, 51.61041647955], + [0.0020847733500204, 51.61041647955], + [-0.0062485600499826, 51.61041647955], ) ) - # expect a RuntimeError: Expected even number of crossings on gridline. snail.core.intersections.split_polygon( bad_poly, 36082,