[geos-commits] [SCM] GEOS branch main updated. 06da4f3a7808288b3aef036aed02c866038985ad
git at osgeo.org
git at osgeo.org
Sun Nov 8 09:16:36 PST 2026
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 06da4f3a7808288b3aef036aed02c866038985ad (commit)
from 2c8e04faf8235456b417979bcf23271ef34148a7 (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 06da4f3a7808288b3aef036aed02c866038985ad
Author: Alejandro <ale at alepacheco.dev>
Date: Mon Nov 9 00:16:01 2026 +0700
GridIntersection: keep traversal on the grid when a segment's slope underflows (#1593)
* GridIntersection: keep traversal on the grid when a segment's slope underflows
crossing() picked the side through which a segment leaves a cell using
the segment's slope dy/dx. For a nearly horizontal segment the slope
underflows to zero (e.g. 2e-30 / 1e300), and for a nearly vertical one
it overflows to infinity, so the wrong side can be chosen. In a padding
cell this sends traverseCells() outside the cell matrix, where
get_cell() only asserts.
- crossing(): interpolate using the fraction of the segment's extent in
x or y instead of its slope. The coordinate is still computed from
the same endpoint, so a segment gives the same result in either
orientation.
- get_cell(): throw instead of asserting when the traversal leaves the
grid.
* Edit comments
---------
Co-authored-by: Dan Baston <dbaston at gmail.com>
diff --git a/src/operation/grid/Cell.cpp b/src/operation/grid/Cell.cpp
index 45b652d7c..f48f728ae 100644
--- a/src/operation/grid/Cell.cpp
+++ b/src/operation/grid/Cell.cpp
@@ -64,7 +64,10 @@ crossing(const geom::Envelope& e, const CoordinateXY& c1, const CoordinateXY& c2
}
}
- const double m = std::abs((c2.y - c1.y) / (c2.x - c1.x));
+ // Interpolate with the fraction of the segment's extent in x or y rather than with its slope, which
+ // underflows to zero or overflows to infinity when the segment is nearly horizontal or vertical.
+ const double dx = std::abs(c2.x - c1.x);
+ const double dy = std::abs(c2.y - c1.y);
const bool up = c2.y > c1.y;
const bool right = c2.x > c1.x;
@@ -73,22 +76,22 @@ crossing(const geom::Envelope& e, const CoordinateXY& c1, const CoordinateXY& c2
if (right) {
// 1st quadrant
- const double y2 = c1.y + m * (e.getMaxX() - c1.x);
+ const double y2 = c1.y + dy * ((e.getMaxX() - c1.x) / dx);
if (y2 < e.getMaxY()) {
return Crossing{ Side::RIGHT, e.getMaxX(), std::clamp(y2, e.getMinY(), e.getMaxY()) };
} else {
- double x2 = c1.x + (e.getMaxY() - c1.y) / m;
+ double x2 = c1.x + dx * ((e.getMaxY() - c1.y) / dy);
return Crossing{ Side::TOP, std::clamp(x2, e.getMinX(), e.getMaxX()), e.getMaxY() };
}
} else {
// 2nd quadrant
- const double y2 = c1.y + m * (c1.x - e.getMinX());
+ const double y2 = c1.y + dy * ((c1.x - e.getMinX()) / dx);
if (y2 < e.getMaxY()) {
return Crossing{ Side::LEFT, e.getMinX(), std::clamp(y2, e.getMinY(), e.getMaxY()) };
} else {
- double x2 = c1.x - (e.getMaxY() - c1.y) / m;
+ double x2 = c1.x - dx * ((e.getMaxY() - c1.y) / dy);
return Crossing{ Side::TOP, std::clamp(x2, e.getMinX(), e.getMaxX()), e.getMaxY() };
}
}
@@ -100,22 +103,22 @@ crossing(const geom::Envelope& e, const CoordinateXY& c1, const CoordinateXY& c2
if (right) {
// 4th quadrant
- const double y2 = c2.y + m * (c2.x - e.getMaxX());
+ const double y2 = c2.y + dy * ((c2.x - e.getMaxX()) / dx);
if (y2 > e.getMinY()) {
return Crossing{ Side::RIGHT, e.getMaxX(), std::clamp(y2, e.getMinY(), e.getMaxY()) };
} else {
- double x2 = c2.x - (e.getMinY() - c2.y) / m;
+ double x2 = c2.x - dx * ((e.getMinY() - c2.y) / dy);
return Crossing{ Side::BOTTOM, std::clamp(x2, e.getMinX(), e.getMaxX()), e.getMinY() };
}
} else {
// 3rd quadrant
- const double y2 = c2.y + m * (e.getMinX() - c2.x);
+ const double y2 = c2.y + dy * ((e.getMinX() - c2.x) / dx);
if (y2 > e.getMinY()) {
return Crossing{ Side::LEFT, e.getMinX(), std::clamp(y2, e.getMinY(), e.getMaxY()) };
} else {
- double x2 = c2.x + (e.getMinY() - c2.y) / m;
+ double x2 = c2.x + dx * ((e.getMinY() - c2.y) / dy);
return Crossing{ Side::BOTTOM, std::clamp(x2, e.getMinX(), e.getMaxX()), e.getMinY() };
}
}
diff --git a/src/operation/grid/GridIntersection.cpp b/src/operation/grid/GridIntersection.cpp
index ef94ce5ad..8440b84fd 100644
--- a/src/operation/grid/GridIntersection.cpp
+++ b/src/operation/grid/GridIntersection.cpp
@@ -69,8 +69,12 @@ GridIntersection::getIntersectionFractions(const Grid<bounded_extent>& raster_gr
static Cell*
get_cell(Matrix<std::unique_ptr<Cell>>& cells, const Grid<infinite_extent>& ex, size_t row, size_t col)
{
- assert(row < cells.getNumRows());
- assert(col < cells.getNumCols());
+ // A traversal that leaves the grid through its first row or column
+ // wraps around to SIZE_MAX, so one comparison covers both directions.
+ if (row >= cells.getNumRows() || col >= cells.getNumCols()) {
+ throw util::GEOSException("Grid traversal left the grid.");
+ }
+
if (cells(row, col) == nullptr) {
cells(row, col) = std::make_unique<Cell>(ex.getCellEnvelope(row, col));
}
diff --git a/tests/unit/operation/grid/GridIntersectionTest.cpp b/tests/unit/operation/grid/GridIntersectionTest.cpp
index 940d54f3a..d1ce00c99 100644
--- a/tests/unit/operation/grid/GridIntersectionTest.cpp
+++ b/tests/unit/operation/grid/GridIntersectionTest.cpp
@@ -1092,4 +1092,66 @@ void object::test<58>()
}
}
+template<>
+template<>
+void object::test<59>()
+{
+ set_test_name("subdividePolygon: traversal stays in grid when the slope of a segment underflows");
+
+ // The first edge crosses y = 0 far to the right of the grid, so it leaves the right
+ // padding cell through its top. Its slope (2e-30 / 1e300) underflows to zero.
+ Envelope e(0, 10, -5, 5);
+ Grid<bounded_extent> ext(e, 1, 5);
+
+ auto g = wkt_reader_.read("POLYGON ((2 -1e-30, 1e300 1e-30, 1e300 1, 2 1, 2 -1e-30))");
+
+ auto subd = GridIntersection::subdividePolygon(ext, *g, false);
+
+ // Part of the polygon inside the grid extent
+ // input manually because g->getFactory()->toGeometry(&e)->intersection(g.get()) gives POLYGON EMPTY
+ auto inside = wkt_reader_.read("POLYGON ((2 0, 10 0, 10 1, 2 1, 2 0))");
+ check_subdivided_polygon(*inside, *subd);
+}
+
+template<>
+template<>
+void object::test<60>()
+{
+ set_test_name("getIntersectionFractions: traversal stays in grid when the slope of a segment underflows");
+
+ // Same as the previous test, with an extra vertex so that the ring is not processed as a rectangle.
+ Envelope e(0, 10, -5, 5);
+ Grid<bounded_extent> ext(e, 1, 5);
+
+ auto g = wkt_reader_.read("POLYGON ((2 -1e-30, 1e300 1e-30, 1e300 1, 5 1, 2 1, 2 -1e-30))");
+
+ auto fracs = GridIntersection::getIntersectionFractions(ext, *g);
+
+ // Part of the polygon inside the grid extent
+ // input manually because g->getFactory()->toGeometry(&e)->intersection(g.get()) gives POLYGON EMPTY
+ auto inside = wkt_reader_.read("POLYGON ((2 0, 10 0, 10 1, 2 1, 2 0))");
+ check_area(*fracs, ext, *inside);
+}
+
+template<>
+template<>
+void object::test<61>()
+{
+ set_test_name("subdividePolygon: traversal stays in grid when the slope of a segment overflows");
+
+ // The first edge crosses x = 0 far below the grid, so it leaves the bottom padding
+ // cell through its left side. Its slope (1e300 / 2e-30) overflows to infinity.
+ Envelope e(-5, 5, -10, 0);
+ Grid<bounded_extent> ext(e, 5, 1);
+
+ auto g = wkt_reader_.read("POLYGON ((1e-30 -2.5, -1e-30 -1e300, -1 -1e300, -1 -2.5, 1e-30 -2.5))");
+
+ auto subd = GridIntersection::subdividePolygon(ext, *g, false);
+
+ // Part of the polygon inside the grid extent
+ // input manually because g->getFactory()->toGeometry(&e)->intersection(g.get()) gives a LINESTRING
+ auto inside = wkt_reader_.read("POLYGON ((-1 -10, 0 -10, 0 -2.5, -1 -2.5, -1 -10))");
+ check_subdivided_polygon(*inside, *subd);
+}
+
}
-----------------------------------------------------------------------
Summary of changes:
src/operation/grid/Cell.cpp | 21 ++++----
src/operation/grid/GridIntersection.cpp | 8 ++-
tests/unit/operation/grid/GridIntersectionTest.cpp | 62 ++++++++++++++++++++++
3 files changed, 80 insertions(+), 11 deletions(-)
hooks/post-receive
--
GEOS
More information about the geos-commits
mailing list