Add grid subdivision / rasterization algorithms

capi/geos_c.cpp

@@ -727,6 +727,13 @@ extern "C" {
         return GEOSClipByRect_r(handle, g, xmin, ymin, xmax, ymax);
     }
 
+    Geometry*
+    GEOSSubdivideByGrid(const Geometry* g, double xmin, double ymin, double xmax, double ymax,
+                        unsigned nx, unsigned ny, int include_exterior)
+    {
+        return GEOSSubdivideByGrid_r(handle, g, xmin, ymin, xmax, ymax, nx, ny, include_exterior);
+    }
+
     int
     GEOSGridIntersectionFractions(const Geometry* g, double xmin, double ymin, double xmax, double ymax,
                                   unsigned nx, unsigned ny, float* buf)

capi/geos_c.h.in

@@ -1168,6 +1168,15 @@ extern GEOSGeometry GEOS_DLL *GEOSClipByRect_r(
     double xmin, double ymin,
     double xmax, double ymax);
 
+/** \see GEOSSubdivideByGrid */
+extern GEOSGeometry GEOS_DLL *GEOSSubdivideByGrid_r(
+    GEOSContextHandle_t handle,
+    const GEOSGeometry* g,
+    double xmin, double ymin,
+    double xmax, double ymax,
+    unsigned nx, unsigned ny,
+    int includeExterior);
+
 /** \see GEOSGridIntersectionFractions */
 extern int GEOS_DLL GEOSGridIntersectionFractions_r(
     GEOSContextHandle_t handle,
@@ -4021,6 +4030,36 @@ extern GEOSGeometry GEOS_DLL *GEOSClipByRect(
     double xmin, double ymin,
     double xmax, double ymax);
 
+/**
+* Compute the intersection of a geometry with each polygon in
+* a rectangular grid.
+* \param g The input geometry to be clipped
+* \param xmin Left bound of grd
+* \param ymin Lower bound of grid
+* \param xmax Right bound of grid
+* \param ymax Upper bound of grid
+* \param nx number of columns in grid
+* \param ny number of rows in grid
+* \param includeExterior whether to include portions of the
+         input geometry that do not intersect the grid in
+         the returned result.
+* \return A geometry collection whose components represent the
+*         intersection with each cell in the grid.
+* Caller is responsible for freeing with GEOSGeom_destroy().
+*
+* This algorithm operates in 2D only; Z and M values in the input will be dropped.
+
+* \see GEOSGetGridIntersectionFractions
+*
+* \since 3.14.0
+*/
+extern GEOSGeometry GEOS_DLL *GEOSSubdivideByGrid(
+    const GEOSGeometry* g,
+    double xmin, double ymin,
+    double xmax, double ymax,
+    unsigned nx, unsigned ny,
+    int includeExterior);
+
 /**
 * Determine the fraction of each cell in a rectangular grid
 * that is covered by a polygon.
@@ -4069,7 +4108,6 @@ extern GEOSGeometry GEOS_DLL *GEOSSharedPaths(
 
 /* ========== Clustering functions ========== */
 /** @name Clustering
-* Functions for clustering geometries
 */
 
 

capi/geos_ts_c.cpp

@@ -1819,6 +1819,20 @@ extern "C" {
         });
     }
 
+    Geometry*
+    GEOSSubdivideByGrid_r(GEOSContextHandle_t extHandle, const Geometry* g, double xmin, double ymin,
+                          double xmax, double ymax, unsigned nx, unsigned ny, int include_exterior)
+    {
+        return execute(extHandle, [&]() {
+            Envelope env(xmin, xmax, ymin, ymax);
+            double dx = env.getWidth() / static_cast<double>(nx);
+            double dy = env.getHeight() / static_cast<double>(ny);
+            geos::operation::grid::Grid<geos::operation::grid::bounded_extent> grid(env, dx, dy);
+
+            return geos::operation::grid::GridIntersection::subdividePolygon(grid, *g, include_exterior).release();
+        });
+    }
+
     int
     GEOSGridIntersectionFractions_r(GEOSContextHandle_t extHandle, const Geometry* g, double xmin, double ymin,
                                     double xmax, double ymax, unsigned nx, unsigned ny, float* buf)

include/geos/geom/Envelope.h

@@ -562,6 +562,19 @@ class GEOS_DLL Envelope {
         return covers(x, y);
     }
 
+    /** \brief
+     * Check if the point p lies in the interior of this Envelope.
+     *
+     * @param p the Coordinate to be tested
+     * @return true if the point lies in the interior
+     */
+    bool
+    containsProperly(const CoordinateXY& p) const
+    {
+        return std::isless(p.x, maxx) && std::isgreater(p.x, minx) &&
+                std::isless(p.y, maxy) && std::isgreater(p.y,  miny);
+    }
+
     /** \brief
      * Check if the point p intersects (lies inside) the region of this Envelope.
      *

include/geos/operation/grid/Cell.h

@@ -40,6 +40,9 @@ class Cell
     {
     }
 
+    /// Get all points that fall on the specified side of this cell.
+    void getEdgePoints(Side s, std::vector<geom::CoordinateXY>& points) const;
+
     const geom::Envelope& box() const { return m_box; }
 
     double getWidth() const;
@@ -76,13 +79,18 @@ class Cell
      * @param c             Coordinate to process
      * @param prev_original The last *uninterpolated* coordinate preceding `c` in the
      *                      boundary being processed
+     * @param exitOnBoundary Whether to exit the cell if `c` lies on the boundary. This is needed to ensure that
+     *                       exit nodes are added to both polygons whose shared boundary is coincident with the
+     *                       Cell boundary. If these nodes are not needed (for example, in the area-only case),
+     *                       it is more efficient to _not_ exit the cell.
+     * @param parentage an optional pointer indicating the source of the coordinate.
      *
      * @return `true` if the Coordinate was inside this cell, `false` otherwise
      */
-    bool take(const geom::CoordinateXY& c, const geom::CoordinateXY* prev_original = nullptr);
+    bool take(const geom::CoordinateXY& c, const geom::CoordinateXY* prev_original, bool exitOnBoundary, const void* parentage);
 
   private:
-    std::vector<const std::vector<geom::CoordinateXY>*> getCoordLists() const;
+    std::vector<const Traversal*> getTraversals() const;
 
     enum class Location
     {
@@ -95,9 +103,9 @@ class Cell
 
     std::vector<Traversal> m_traversals;
 
-    Side side(const geom::CoordinateXY& c) const;
+    Side getSide(const geom::CoordinateXY& c) const;
 
-    Location location(const geom::CoordinateXY& c) const;
+    Location getLocation(const geom::CoordinateXY& c) const;
 
     /// If no Traversals have been started or the most recent Traversal has been completed,
     /// return a new Traversal. Otherwise, return the most recent Traversal.

include/geos/operation/grid/Grid.h

@@ -20,8 +20,6 @@
 
 #include <geos/geom/Envelope.h>
 
-constexpr double DEFAULT_GRID_COMPAT_TOL = 1e-6;
-
 namespace geos::operation::grid {
 struct infinite_extent
 {
@@ -33,26 +31,40 @@ struct bounded_extent
     static constexpr size_t padding = 0;
 };
 
+/**
+ * @brief The Grid class represents a grid of constant-size rectangular cells that covers a specified envelope.
+ *        The width of the cells may be different from the height. If the Grid has an "infinite" extent rather than
+ *        a "bounded" extent, then an extra row and column will be added on all side of the grid. The size of the cells
+ *        in these columns may be larger than those in the primary grid, such that the extended grid covers a "domain"
+ *        that is larger than the extent of the regular grid.
+ */
 template<typename extent_tag>
-class Grid
+class GEOS_DLL Grid
 {
 
   public:
+    /// Construct a bounded grid covering a specified extent.
     Grid(const geom::Envelope& extent, double dx, double dy)
       : m_extent{ extent }
       , m_domain{}
       , m_dx{ dx }
       , m_dy{ dy }
+      , m_xOrig { extent.getMinX() }
+      , m_yOrig { extent.getMaxY() }
       , m_num_rows{ 2 * extent_tag::padding + (extent.getMaxY() > extent.getMinY() ? static_cast<size_t>(std::round(extent.getHeight() / dy)) : 0) }
       , m_num_cols{ 2 * extent_tag::padding + (extent.getMaxX() > extent.getMinX() ? static_cast<size_t>(std::round(extent.getWidth() / dx)) : 0) } {
         static_assert(std::is_same_v<extent_tag, bounded_extent>);
     }
 
+    /// Construct an infinite grid covering a specified extent with regularly-sized cells, and adding a row and column
+    /// of variably-sized cells to each edge of the primary grid such that the specified domain is covered.
     Grid(const geom::Envelope& extent, double dx, double dy, const geom::Envelope& domain)
       : m_extent{ extent }
       , m_domain{ domain }
       , m_dx{ dx }
       , m_dy{ dy }
+      , m_xOrig { extent.getMinX() }
+      , m_yOrig { extent.getMaxY() }
       , m_num_rows{ 2 * extent_tag::padding + (extent.getMaxY() > extent.getMinY() ? static_cast<size_t>(std::round(extent.getHeight() / dy)) : 0) }
       , m_num_cols{ 2 * extent_tag::padding + (extent.getMaxX() > extent.getMinX() ? static_cast<size_t>(std::round(extent.getWidth() / dx)) : 0) }
     {
@@ -69,6 +81,7 @@ class Grid
         }
     }
 
+    /// Get the column in which the specified x coordinate would fall.
     size_t getColumn(double x) const
     {
         if (extent_tag::padding) {
@@ -94,6 +107,7 @@ class Grid
           getColumn(m_extent.getMaxX()));
     }
 
+    /// Get the row in which the specified y coordinate would fall.
     size_t getRow(double y) const
     {
         if (extent_tag::padding) {
@@ -119,6 +133,8 @@ class Grid
           getRow(m_extent.getMinY()));
     }
 
+    /// Get the cell index in which the specified x and y values would fall. Cells are indexed from left-to-right,
+    /// then top-to-bottom.
     std::size_t getCell(double x, double y) const
     {
         return getRow(y) * getNumCols() + getColumn(x);
@@ -146,15 +162,36 @@ class Grid
 
     const geom::Envelope& getExtent() const { return m_extent; }
 
+    /// Return the number of rows by which another grid is offset from this Grid. It is assumed that the two
+    /// grids have the same resolution, and that the maximum y value of the other grid is less than or equal
+    /// to the maximum y value of this grid.
     size_t getRowOffset(const Grid& other) const { return static_cast<size_t>(std::round(std::abs(other.m_extent.getMaxY() - m_extent.getMaxY()) / m_dy)); }
 
+    /// Return the number of columns by which another grid is offset from this Grid. It is assumed that the two
+    /// grids have the same resolution, and that the minimum x value of the other grid is greater than or equal
+    /// to the minimum x value of this grid.
     size_t getColOffset(const Grid& other) const { return static_cast<size_t>(std::round(std::abs(m_extent.getMinX() - other.m_extent.getMinX()) / m_dx)); }
 
+    /// Get the x coordinate at the center of the specified column.
     double getColX(size_t col) const { return m_extent.getMinX() + (static_cast<double>(col - extent_tag::padding) + 0.5) * m_dx; }
 
+    /// Get the y coordinate at the center of the specified row.
     double getRowY(size_t row) const { return m_extent.getMaxY() - (static_cast<double>(row - extent_tag::padding) + 0.5) * m_dy; }
 
-    Grid<extent_tag> shrinkToFit(const geom::Envelope& b) const
+    Grid crop(const geom::Envelope& e) const
+    {
+        if (m_extent.intersects(e)) {
+            return shrinkToFit(e.intersection(m_extent));
+        } else {
+            return makeEmpty();
+        }
+    }
+
+    /// Reduce the size of the grid to contain only the provided Envelope
+    /// If calcExtentFromNewGrid is true, then the xmax and ymin of the new
+    /// grid will be calculated relative to the origin point of the original grid
+    /// rather than the newly cropped grid.
+    Grid<extent_tag> shrinkToFit(const geom::Envelope& b, bool calcExtentFromNewGrid=true) const
     {
         if (b.getArea() == 0) {
             return makeEmpty();
@@ -171,8 +208,8 @@ class Grid
         double snapped_xmin = m_extent.getMinX() + static_cast<double>(col0 - extent_tag::padding) * m_dx;
         double snapped_ymax = m_extent.getMaxY() - static_cast<double>(row1 - extent_tag::padding) * m_dy;
 
-        // Make sure x0 and y1 are within the reduced extent. Because of
-        // floating point round-off errors, this is not always the case.
+        // Make sure snapped_xmin and snapped_ymax are within the reduced extent.
+        // Because of floating point round-off errors, this is not always the case.
         if (b.getMinX() < snapped_xmin) {
             snapped_xmin -= m_dx;
             col0--;
@@ -201,14 +238,17 @@ class Grid
             num_cols--;
         }
 
+        const double snapped_xmax = m_extent.getMinX() + static_cast<double>(col0 + num_cols - extent_tag::padding) * m_dx;
+        const double snapped_ymin = m_extent.getMaxY() - static_cast<double>(row1 + num_rows - extent_tag::padding) * m_dy;
+
         // Perform offsets relative to the new xmin, ymax origin
         // points. If this is not done, then floating point roundoff
         // error can cause progressive shrink() calls with the same
         // inputs to produce different results.
         geom::Envelope reduced_box = {
             snapped_xmin,
-            std::max(snapped_xmin + static_cast<double>(num_cols) * m_dx, b.getMaxX()),
-            std::min(snapped_ymax - static_cast<double>(num_rows) * m_dy, b.getMinY()),
+            calcExtentFromNewGrid ? std::max(snapped_xmin + static_cast<double>(num_cols) * m_dx, b.getMaxX()) : snapped_xmax,
+            calcExtentFromNewGrid ? std::min(snapped_ymax - static_cast<double>(num_rows) * m_dy, b.getMinY()) : snapped_ymin,
             snapped_ymax
         };
 
@@ -232,22 +272,29 @@ class Grid
         }
 
         if constexpr (std::is_same_v<extent_tag, bounded_extent>) {
-            Grid<extent_tag> reduced{ reduced_box, m_dx, m_dy };
-
-        if (b.getMinX() < reduced.xmin() || b.getMinY() < reduced.ymin() || b.getMaxX() > reduced.xmax() || b.getMaxY() > reduced.ymax()) {
-            throw std::runtime_error("Shrink operation failed.");
-        }
-
+            Grid<extent_tag> reduced{reduced_box, m_dx, m_dy};
+
+            if (!calcExtentFromNewGrid) {
+                reduced.m_xOrig = m_xOrig;
+                reduced.m_yOrig = m_yOrig;
+                reduced.m_skipRows = reduced.getRowOffset(*this);
+                reduced.m_skipCols = reduced.getColOffset(*this);
+            } else if (!reduced.getExtent().contains(b)) {
+                throw std::runtime_error("Shrink operation failed.");
+            }
             return reduced;
         } else {
-            Grid<extent_tag> reduced{ reduced_box, m_dx, m_dy , m_domain };
-
-        if (b.getMinX() < reduced.xmin() || b.getMinY() < reduced.ymin() || b.getMaxX() > reduced.xmax() || b.getMaxY() > reduced.ymax()) {
-            throw std::runtime_error("Shrink operation failed.");
-        }
-
+            Grid<extent_tag> reduced{reduced_box, m_dx, m_dy, m_domain};
+
+            if (!calcExtentFromNewGrid) {
+                reduced.m_xOrig = m_xOrig;
+                reduced.m_yOrig = m_yOrig;
+                reduced.m_skipRows = reduced.getRowOffset(*this);
+                reduced.m_skipCols = reduced.getColOffset(*this);
+            }else if (!reduced.getExtent().contains(b)) {
+                throw std::runtime_error("Shrink operation failed.");
+            }
             return reduced;
-
         }
     }
 
@@ -272,8 +319,18 @@ class Grid
     double m_dx;
     double m_dy;
 
+    double m_xOrig; // origin point that is distinct from xmin of m_extent. Used to allow a subgrid to calculate
+                    // sub-envelopes that exactly match those of the parent grid
+    double m_yOrig;
+
+    size_t m_skipRows{0}; // number of rows to skip when computing a cell envelope using m_yOrig
+    size_t m_skipCols{0}; // number of cols to skip when computing a cell envelope using m_xOrig
+
     size_t m_num_rows;
     size_t m_num_cols;
+
+    friend Grid<infinite_extent> make_infinite(const Grid<bounded_extent>&, const geom::Envelope&);
+    friend Grid<bounded_extent> make_finite(const Grid<infinite_extent>&);
 };
 
 Grid<infinite_extent>