From fd501cb60577d1bfe65faf6081dfbdcdd9e5b26b Mon Sep 17 00:00:00 2001
From: Tom Russell <tom.russell@ouce.ox.ac.uk>
Date: Fri, 24 May 2024 17:56:32 +0100
Subject: [PATCH] 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<linestr> splitAlongGridlines(linestr exterior_crossings,
 
   std::vector<geometry::Coord> crossings_on_gridline;
   std::vector<linestr> 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<linestr> 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<linestr> 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 <class T>
 typename std::enable_if<!std::numeric_limits<T>::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<T>::epsilon() * std::fabs(x + y) * ulp
-         // unless the result is subnormal
-         || std::fabs(x - y) < std::numeric_limits<T>::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<T>::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<T>::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,