Merge pull request #62 from nismod/fix/linestring

Fix linestring splitting issues

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 {

extension/src/grid.hpp

@@ -29,14 +29,14 @@ struct Grid {
   std::vector<double> data;
 
   Grid()
-      : ncols{0}, nrows{0},
-        grid_to_world{Affine()}, data{std::vector<double>()} {
+      : ncols{0}, nrows{0}, grid_to_world{Affine()},
+        data{std::vector<double>()} {
     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<double>()} {
+      : ncols{ncols}, nrows{nrows}, grid_to_world{grid_to_world},
+        data{std::vector<double>()} {
     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);
@@ -155,7 +153,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;

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;
     }
 

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

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",
     [
@@ -58,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,