diff --git a/capi/geos_c.cpp b/capi/geos_c.cpp index 35d14fd18..d26db0966 100644 --- a/capi/geos_c.cpp +++ b/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) diff --git a/capi/geos_c.h.in b/capi/geos_c.h.in index 7c80f7016..942663671 100644 --- a/capi/geos_c.h.in +++ b/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 */ diff --git a/capi/geos_ts_c.cpp b/capi/geos_ts_c.cpp index f4690f04f..d8b57fd87 100644 --- a/capi/geos_ts_c.cpp +++ b/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(nx); + double dy = env.getHeight() / static_cast(ny); + geos::operation::grid::Grid 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) diff --git a/include/geos/geom/Envelope.h b/include/geos/geom/Envelope.h index 61e074be8..687dbbae1 100644 --- a/include/geos/geom/Envelope.h +++ b/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. * diff --git a/include/geos/operation/grid/Cell.h b/include/geos/operation/grid/Cell.h index a58ffc8f0..a307cdaf4 100644 --- a/include/geos/operation/grid/Cell.h +++ b/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& 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*> getCoordLists() const; + std::vector getTraversals() const; enum class Location { @@ -95,9 +103,9 @@ class Cell std::vector 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. diff --git a/include/geos/operation/grid/Grid.h b/include/geos/operation/grid/Grid.h index b18b84b5c..d94d80916 100644 --- a/include/geos/operation/grid/Grid.h +++ b/include/geos/operation/grid/Grid.h @@ -20,8 +20,6 @@ #include -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 -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(std::round(extent.getHeight() / dy)) : 0) } , m_num_cols{ 2 * extent_tag::padding + (extent.getMaxX() > extent.getMinX() ? static_cast(std::round(extent.getWidth() / dx)) : 0) } { static_assert(std::is_same_v); } + /// 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(std::round(extent.getHeight() / dy)) : 0) } , m_num_cols{ 2 * extent_tag::padding + (extent.getMaxX() > extent.getMinX() ? static_cast(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(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(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(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(row - extent_tag::padding) + 0.5) * m_dy; } - Grid 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 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(col0 - extent_tag::padding) * m_dx; double snapped_ymax = m_extent.getMaxY() - static_cast(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(col0 + num_cols - extent_tag::padding) * m_dx; + const double snapped_ymin = m_extent.getMaxY() - static_cast(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(num_cols) * m_dx, b.getMaxX()), - std::min(snapped_ymax - static_cast(num_rows) * m_dy, b.getMinY()), + calcExtentFromNewGrid ? std::max(snapped_xmin + static_cast(num_cols) * m_dx, b.getMaxX()) : snapped_xmax, + calcExtentFromNewGrid ? std::min(snapped_ymax - static_cast(num_rows) * m_dy, b.getMinY()) : snapped_ymin, snapped_ymax }; @@ -232,22 +272,29 @@ class Grid } if constexpr (std::is_same_v) { - Grid 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 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 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 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 make_infinite(const Grid&, const geom::Envelope&); + friend Grid make_finite(const Grid&); }; Grid diff --git a/include/geos/operation/grid/GridIntersection.h b/include/geos/operation/grid/GridIntersection.h index e85d401cb..c97d4e8b5 100644 --- a/include/geos/operation/grid/GridIntersection.h +++ b/include/geos/operation/grid/GridIntersection.h @@ -45,13 +45,13 @@ class GEOS_DLL GridIntersection * @brief Compute the fraction of each cell in a rectangular grid that is covered by a Geometry. * A matrix can be provided to which the fractions will be added. */ - GridIntersection(const Grid& raster_grid, const geom::Geometry& g, std::shared_ptr> cov = nullptr); + GridIntersection(const Grid& raster_grid, const geom::Geometry& g, const std::shared_ptr>& cov = nullptr); /** * @brief Compute the fraction of each cell in a rectangular grid that is covered by an Envelope. * A matrix can be provided to which the fractions will be added. */ - GridIntersection(const Grid& raster_grid, const geom::Envelope& box, std::shared_ptr> cov = nullptr); + GridIntersection(const Grid& raster_grid, const geom::Envelope& box, const std::shared_ptr>& cov = nullptr); /** * @brief Return the intersection result matrix @@ -60,14 +60,25 @@ class GEOS_DLL GridIntersection /** * @brief Partition a polygonal geometry by a grid + * @param grid the Grid by which the polygon should be partitioned. The Grid does not need to cover the polygon. + * @param g a polygonal geometry + * @param includeExterior whether the result should contain any portions of the input that fall outside the grid. + * @return a GeometryCollection containing the partitioned polygon. */ - static std::unique_ptr subdividePolygon(const Grid& p_grid, const geom::Geometry& g, bool includeExterior); + static std::unique_ptr subdividePolygon(const Grid& grid, const geom::Geometry& g, bool includeExterior); + /** + * @brief Calculate the fraction of each cell in a Grid that is covered by a polygon. + * @param grid the Grid for which coverage fractions should be returned + * @param g a polygonal or linear geometry + * @return a matrix having the same number of rows and columns as the grid, whose values contain either the + * fraction of the cell area (for polygonal inputs) or the length of the intersection (for linear inputs) + */ static std::shared_ptr> - getIntersectionFractions(const Grid& raster_grid, const geom::Geometry& g); + getIntersectionFractions(const Grid& grid, const geom::Geometry& g); static std::shared_ptr> - getIntersectionFractions(const Grid& raster_grid, const geom::Envelope& box); + getIntersectionFractions(const Grid& grid, const geom::Envelope& box); /** * @brief Determines the bounding box of the raster-vector intersection. Considers the bounding boxes diff --git a/include/geos/operation/grid/PerimeterDistance.h b/include/geos/operation/grid/PerimeterDistance.h index 86691d257..b02ac2750 100644 --- a/include/geos/operation/grid/PerimeterDistance.h +++ b/include/geos/operation/grid/PerimeterDistance.h @@ -38,12 +38,23 @@ class GEOS_DLL PerimeterDistance { /** Calculates the counter-clockwise distance between two locations on the perimeter * of an Envelope. The locations are specified as clockwise distances from the - * lower left corner of the Envelope, consistent with the `perimeter_distance` + * lower left corner of the Envelope, consistent with the `getPerimeterDistance` * function. */ static double getPerimeterDistanceCCW(double measure1, double measure2, double perimeter); -}; + /** Tests whether measure x would be passed when traveling from a to b in a CCW + * direction. + */ + static bool + isBetweenCCW(double x, double a, double b); + /** Tests whether the perimeter distance of c1 is less than the perimeter distance of c2, in a + * robust way. + */ + static bool + isLessThan(const geom::Envelope& e, const geom::CoordinateXY& c1, const geom::CoordinateXY& c2); + +}; } diff --git a/include/geos/operation/grid/Side.h b/include/geos/operation/grid/Side.h index b9015f7a1..f91866908 100644 --- a/include/geos/operation/grid/Side.h +++ b/include/geos/operation/grid/Side.h @@ -27,7 +27,28 @@ enum class Side BOTTOM }; -std::ostream& -operator<<(std::ostream& os, const Side& s); +inline std::ostream& +operator<<(std::ostream& os, const Side& s) +{ + switch (s) { + case Side::NONE: + os << "none"; + return os; + case Side::LEFT: + os << "left"; + return os; + case Side::RIGHT: + os << "right"; + return os; + case Side::TOP: + os << "top"; + return os; + case Side::BOTTOM: + os << "bottom"; + return os; + } + + return os; // unreachable statement needed for -Werror=return-type +} } diff --git a/include/geos/operation/grid/Traversal.h b/include/geos/operation/grid/Traversal.h index 68bc516d7..67164cc0b 100644 --- a/include/geos/operation/grid/Traversal.h +++ b/include/geos/operation/grid/Traversal.h @@ -14,11 +14,12 @@ #pragma once -#include - +#include #include #include +#include + namespace geos::operation::grid { /** @@ -26,17 +27,20 @@ namespace geos::operation::grid { * within a grid cell, as well as the `Side` from which the line * entered and exited the cell. */ -class Traversal +class GEOS_DLL Traversal { public: Traversal() : m_entry{ Side::NONE } , m_exit{ Side::NONE } + , m_parentage{ nullptr } { } bool isClosedRing() const; + bool isClosedRingWithArea() const; + bool isEmpty() const; bool isEntered() const; @@ -48,7 +52,7 @@ class Traversal bool hasMultipleUniqueCoordinates() const; /// Begin a Traversal on the specified `Side` - void enter(const geom::CoordinateXY& c, Side s); + void enter(const geom::CoordinateXY& c, Side s, const void* parentage); /// Complete a Traversal on the specified `Side` void exit(const geom::CoordinateXY& c, Side s); @@ -57,6 +61,8 @@ class Traversal Side getExitSide() const { return m_exit; } + const geom::CoordinateXY& getFirstCoordinate() const; + const geom::CoordinateXY& getLastCoordinate() const; const geom::CoordinateXY& getExitCoordinate() const; @@ -67,10 +73,13 @@ class Traversal const std::vector& getCoordinates() const { return m_coords; } + const void* getParentage() const { return m_parentage; } + private: std::vector m_coords; Side m_entry; Side m_exit; + const void* m_parentage; }; } diff --git a/include/geos/operation/grid/TraversalAreas.h b/include/geos/operation/grid/TraversalAreas.h index 5d660c909..d8f88c327 100644 --- a/include/geos/operation/grid/TraversalAreas.h +++ b/include/geos/operation/grid/TraversalAreas.h @@ -19,6 +19,8 @@ namespace geos::operation::grid { +class Traversal; + class GEOS_DLL TraversalAreas { public: /** @@ -32,7 +34,7 @@ class GEOS_DLL TraversalAreas { * @return total area */ static double - getLeftHandArea(const geom::Envelope& box, const std::vector*>& coord_lists); + getLeftHandArea(const geom::Envelope& box, const std::vector& coord_lists); /** * @brief Return an areal geometry representing the closed rings formed by this box and the provided Coordinate sequences @@ -46,7 +48,7 @@ class GEOS_DLL TraversalAreas { * @return a Polygon or MultiPolygon geometry */ static std::unique_ptr - getLeftHandRings(const geom::GeometryFactory& gfact, const geom::Envelope& box, const std::vector*>& coord_lists); + getLeftHandRings(const geom::GeometryFactory& gfact, const geom::Envelope& box, const std::vector& coord_lists); }; diff --git a/src/operation/grid/Cell.cpp b/src/operation/grid/Cell.cpp index 3ae5cb156..45b652d7c 100644 --- a/src/operation/grid/Cell.cpp +++ b/src/operation/grid/Cell.cpp @@ -19,17 +19,18 @@ #include #include +#define DEBUG_CELL 0 +#include + +#include "geos/geom/GeometryFactory.h" + using geos::geom::CoordinateXY; using geos::geom::Geometry; using geos::geom::GeometryFactory; namespace geos::operation::grid { -static bool -strictly_contains(const geom::Envelope& e, const geom::CoordinateXY& c) { - return e.getMinX() < c.x && e.getMaxX() > c.x && e.getMinY() < c.y && e.getMaxY() > c.y; -} - +// Report where we leave an Envelope e when moving from c1 to c2 static Crossing crossing(const geom::Envelope& e, const CoordinateXY& c1, const CoordinateXY& c2) { @@ -39,6 +40,10 @@ crossing(const geom::Envelope& e, const CoordinateXY& c1, const CoordinateXY& c2 return Crossing{ Side::TOP, c1.x, e.getMaxY() }; } else if (c2.y <= e.getMinY()) { return Crossing{ Side::BOTTOM, c1.x, e.getMinY() }; + //} else if (c2.x == e.getMinX()) { + // return Crossing{ Side::LEFT, c2.x, c2.y }; + //} else if (c2.x == e.getMaxX()) { + // return Crossing{ Side::RIGHT, c2.x, c2.y }; } else { throw std::runtime_error("Never get here."); } @@ -50,20 +55,25 @@ crossing(const geom::Envelope& e, const CoordinateXY& c1, const CoordinateXY& c2 return Crossing{ Side::RIGHT, e.getMaxX(), c1.y }; } else if (c2.x <= e.getMinX()) { return Crossing{ Side::LEFT, e.getMinX(), c1.y }; + //} else if (c2.y == e.getMinY()) { + // return Crossing{ Side::BOTTOM, c2.x, c2.y }; + //} else if (c2.y == e.getMaxY()) { + // return Crossing{ Side::TOP, c2.x, c2.y }; } else { throw std::runtime_error("Never get here"); } } - double m = std::abs((c2.y - c1.y) / (c2.x - c1.x)); + const double m = std::abs((c2.y - c1.y) / (c2.x - c1.x)); - bool up = c2.y > c1.y; - bool right = c2.x > c1.x; + const bool up = c2.y > c1.y; + const bool right = c2.x > c1.x; if (up) { + if (right) { // 1st quadrant - double y2 = c1.y + m * (e.getMaxX() - c1.x); + const double y2 = c1.y + m * (e.getMaxX() - c1.x); if (y2 < e.getMaxY()) { return Crossing{ Side::RIGHT, e.getMaxX(), std::clamp(y2, e.getMinY(), e.getMaxY()) }; @@ -73,7 +83,7 @@ crossing(const geom::Envelope& e, const CoordinateXY& c1, const CoordinateXY& c2 } } else { // 2nd quadrant - double y2 = c1.y + m * (c1.x - e.getMinX()); + const double y2 = c1.y + m * (c1.x - e.getMinX()); if (y2 < e.getMaxY()) { return Crossing{ Side::LEFT, e.getMinX(), std::clamp(y2, e.getMinY(), e.getMaxY()) }; @@ -83,24 +93,29 @@ crossing(const geom::Envelope& e, const CoordinateXY& c1, const CoordinateXY& c2 } } } else { + // For downward segments, we calculate constructed coordinates relative to c2, not c1. This is so the + // same coordinate will be calculated regardless of the segment orientation. This is important for maintaining + // valid polygon coverages (the same segment will be processed with opposite orientation along shared + // boundaries) + if (right) { // 4th quadrant - double y2 = c1.y - m * (e.getMaxX() - c1.x); + const double y2 = c2.y + m * (c2.x - e.getMaxX()); if (y2 > e.getMinY()) { return Crossing{ Side::RIGHT, e.getMaxX(), std::clamp(y2, e.getMinY(), e.getMaxY()) }; } else { - double x2 = c1.x + (c1.y - e.getMinY()) / m; + double x2 = c2.x - (e.getMinY() - c2.y) / m; return Crossing{ Side::BOTTOM, std::clamp(x2, e.getMinX(), e.getMaxX()), e.getMinY() }; } } else { // 3rd quadrant - double y2 = c1.y - m * (c1.x - e.getMinX()); + const double y2 = c2.y + m * (e.getMinX() - c2.x); if (y2 > e.getMinY()) { return Crossing{ Side::LEFT, e.getMinX(), std::clamp(y2, e.getMinY(), e.getMaxY()) }; } else { - double x2 = c1.x - (c1.y - e.getMinY()) / m; + double x2 = c2.x + (e.getMinY() - c2.y) / m; return Crossing{ Side::BOTTOM, std::clamp(x2, e.getMinX(), e.getMaxX()), e.getMinY() }; } } @@ -127,7 +142,7 @@ Cell::getArea() const } Side -Cell::side(const CoordinateXY& c) const +Cell::getSide(const CoordinateXY& c) const { if (c.x == m_box.getMinX()) { return Side::LEFT; @@ -151,15 +166,15 @@ Cell::forceExit() const CoordinateXY& last = getLastTraversal().getLastCoordinate(); - if (location(last) == Location::BOUNDARY) { - getLastTraversal().forceExit(side(last)); + if (getLocation(last) == Location::BOUNDARY) { + getLastTraversal().forceExit(getSide(last)); } } Cell::Location -Cell::location(const CoordinateXY& c) const +Cell::getLocation(const CoordinateXY& c) const { - if (strictly_contains(m_box, c)) { + if (m_box.containsProperly(c)) { return Cell::Location::INSIDE; } @@ -187,19 +202,27 @@ Cell::getLastTraversal() } bool -Cell::take(const CoordinateXY& c, const CoordinateXY* prev_original) +Cell::take(const CoordinateXY& c, const CoordinateXY* prev_original, bool exitOnBoundary, const void* parentage) { Traversal& t = traversal_in_progress(); if (t.isEmpty()) { - // std::cout << "Entering " << m_box << " from " << side(c) << " at " << c << std::endl; +#if DEBUG_CELL + std::cout << std::setprecision(17); + std::cout << "Entering " << m_box << " from " << getSide(c) << " at " << c << std::endl; +#endif - t.enter(c, side(c)); + t.enter(c, getSide(c), parentage); return true; } - if (location(c) != Cell::Location::OUTSIDE) { - // std::cout << "Still in " << m_box << " with " << c << std::endl; + const auto loc = getLocation(c); + const bool canTakeCoordinate = exitOnBoundary ? (loc == Location::INSIDE) : (loc != Location::OUTSIDE); + + if (canTakeCoordinate) { +#if DEBUG_CELL + std::cout << "Still in " << m_box << " with " << c << std::endl; +#endif t.add(c); @@ -208,6 +231,17 @@ Cell::take(const CoordinateXY& c, const CoordinateXY* prev_original) } return true; + } else if (loc == Location::BOUNDARY) { + if (c.x == m_box.getMaxX()) { + t.exit(c, Side::RIGHT); + } else if (c.y == m_box.getMaxY()) { + t.exit(c, Side::TOP); + } else if (c.x == m_box.getMinX()) { + t.exit(c, Side::LEFT); + } else { + t.exit(c, Side::BOTTOM); + } + return false; } // We need to calculate the coordinate of the cell exit point using only uninterpolated coordinates. @@ -217,8 +251,10 @@ Cell::take(const CoordinateXY& c, const CoordinateXY* prev_original) Crossing x = prev_original ? crossing(m_box, *prev_original, c) : crossing(m_box, t.getLastCoordinate(), c); t.exit(x.getCoord(), x.getSide()); - // std::cout << "Leaving " << m_box << " from " << x.side() << " at " << x.coord(); - // std::cout << " on the way to " << c << std::endl; +#if DEBUG_CELL + std::cout << "Leaving " << m_box << " from " << x.getSide() << " at " << x.getCoord(); + std::cout << " on the way to " << c << std::endl; +#endif return false; } @@ -235,7 +271,11 @@ bool Cell::isDetermined() const { for (const auto& t : m_traversals) { - if (!t.isTraversed() && !t.isClosedRing()) { + if (t.isClosedRing()) { + if (!t.isClosedRingWithArea()) { + continue; + } + } else if (!t.isTraversed()) { continue; } @@ -247,33 +287,49 @@ Cell::isDetermined() const return false; } -std::vector*> -Cell::getCoordLists() const +std::vector +Cell::getTraversals() const { - std::vector*> coord_lists; - coord_lists.reserve(m_traversals.size()); + std::vector traversals; + traversals.reserve(m_traversals.size()); for (const auto& t : m_traversals) { - if (t.isTraversed() || t.isClosedRing()) { - coord_lists.push_back(&t.getCoordinates()); + if (t.isTraversed() || t.isClosedRingWithArea()) { + traversals.push_back(&t); } } - return coord_lists; + return traversals; } double Cell::getCoveredFraction() const { - auto coord_lists = getCoordLists(); + auto coord_lists = getTraversals(); return TraversalAreas::getLeftHandArea(m_box, coord_lists) / getArea(); } std::unique_ptr Cell::getCoveredPolygons(const GeometryFactory& gfact) const { - auto coord_lists = getCoordLists(); + auto coord_lists = getTraversals(); return TraversalAreas::getLeftHandRings(gfact, m_box, coord_lists); } +void Cell::getEdgePoints(Side s, std::vector &edgePoints) const +{ + for (const Traversal& t : m_traversals) { + const auto& coords = t.getCoordinates(); + + for (const auto& c : coords) { + if ((s == Side::LEFT && c.x == m_box.getMinX()) || + (s == Side::RIGHT && c.x == m_box.getMaxX()) || + (s == Side::BOTTOM && c.y == m_box.getMinY()) || + (s == Side::TOP && c.y == m_box.getMaxY())) { + edgePoints.push_back(c); + } + } + } +} + } diff --git a/src/operation/grid/Grid.cpp b/src/operation/grid/Grid.cpp index 28b84c709..d4db6659b 100644 --- a/src/operation/grid/Grid.cpp +++ b/src/operation/grid/Grid.cpp @@ -19,24 +19,24 @@ using geos::geom::Envelope; namespace geos::operation::grid { template<> -Envelope +Envelope GEOS_DLL Grid::getCellEnvelope(size_t row, size_t col) const { // The ternary clauses below are used to make sure that the cells along - // the right and bottom edges of our grid are slightly larger than m_dx,dy + // the right and bottom edges of our grid are slightly larger than m_dx, m_dy // if needed to make sure that we capture our whole extent. This is necessary // because xmin + nx*m_dx may be less than xmax because of floating point // errors. return { - xmin() + static_cast(col) * dx(), - row == (getNumRows() - 1) ? ymin() : (ymax() - static_cast(row + 1) * dy()), - col == (getNumCols() - 1) ? xmax() : (xmin() + static_cast(col + 1) * dx()), - ymax() - static_cast(row) * dy() + m_xOrig + static_cast(col + m_skipCols) * dx(), + col == (getNumCols() - 1) ? xmax() : (m_xOrig + static_cast(col + m_skipCols + 1) * dx()), + row == (getNumRows() - 1) ? ymin() : (m_yOrig - static_cast(row + m_skipRows + 1) * dy()), + m_yOrig - static_cast(row + m_skipRows) * dy() }; } template<> -Envelope +Envelope GEOS_DLL Grid::getCellEnvelope(size_t row, size_t col) const { double cell_xmin, cell_xmax, cell_ymin, cell_ymax; @@ -44,56 +44,76 @@ Grid::getCellEnvelope(size_t row, size_t col) const constexpr double PADDED_CELL_SIZE = 1.0; if (col == 0) { + // Left-side padding cell_xmin = std::min(xmin() - PADDED_CELL_SIZE, m_domain.getMinX()); } else if (col == getNumCols() - 1) { + // Right-side padding cell_xmin = xmax(); // because rightmost col of regular may have different width from others } else { - cell_xmin = xmin() + static_cast(col - 1) * dx(); + // Internal + cell_xmin = m_xOrig + static_cast(m_skipCols + col - 1) * dx(); } switch (getNumCols() - col) { case 1: + // Right-side padding cell_xmax = std::max(xmax() + PADDED_CELL_SIZE, m_domain.getMaxX()); break; case 2: + // Rightmost internal cell cell_xmax = xmax(); break; default: - cell_xmax = xmin() + static_cast(col) * dx(); + cell_xmax = m_xOrig + static_cast(m_skipCols + col) * dx(); } if (row == 0) { + // Top row padding cell_ymax = std::max(ymax() + PADDED_CELL_SIZE, m_domain.getMaxY()); } else if (row == getNumRows() - 1) { + // Bottom row padding cell_ymax = ymin(); // because bottom row of regular may have different height from others } else { - cell_ymax = ymax() - static_cast(row - 1) * dy(); + cell_ymax = m_yOrig - static_cast(m_skipRows + row - 1) * dy(); } switch (getNumRows() - row) { case 1: + // Bottom row padding cell_ymin = std::min(xmin() - PADDED_CELL_SIZE, m_domain.getMinY()); break; case 2: + // Bottom-most internal row cell_ymin = ymin(); break; default: - cell_ymin = ymax() - static_cast(row) * dy(); + cell_ymin = m_yOrig - static_cast(m_skipRows + row) * dy(); } - return { cell_xmin, cell_xmax, cell_ymin, cell_ymax }; + Envelope box{ cell_xmin, cell_xmax, cell_ymin, cell_ymax }; + return box; } Grid make_infinite(const Grid& grid, const Envelope& domain) { - return { grid.getExtent(), grid.dx(), grid.dy() , domain }; + Grid ret{ grid.getExtent(), grid.dx(), grid.dy() , domain }; + ret.m_xOrig = grid.m_xOrig; + ret.m_yOrig = grid.m_yOrig; + ret.m_skipCols = grid.m_skipCols; + ret.m_skipRows = grid.m_skipRows; + return ret; } Grid make_finite(const Grid& grid) { - return { grid.getExtent(), grid.dx(), grid.dy() }; + Grid ret{ grid.getExtent(), grid.dx(), grid.dy() }; + ret.m_xOrig = grid.m_xOrig; + ret.m_yOrig = grid.m_yOrig; + ret.m_skipCols = grid.m_skipCols; + ret.m_skipRows = grid.m_skipRows; + return ret; } } diff --git a/src/operation/grid/GridIntersection.cpp b/src/operation/grid/GridIntersection.cpp index 54e67424f..7c4815b75 100644 --- a/src/operation/grid/GridIntersection.cpp +++ b/src/operation/grid/GridIntersection.cpp @@ -22,6 +22,7 @@ #include #include #include +#include #include #include @@ -50,6 +51,8 @@ GridIntersection::getIntersectionFractions(const Grid& raster_gr static Cell* get_cell(Matrix>& cells, const Grid& ex, size_t row, size_t col) { + assert(row < cells.getNumRows()); + assert(col < cells.getNumCols()); if (cells(row, col) == nullptr) { cells(row, col) = std::make_unique(ex.getCellEnvelope(row, col)); } @@ -82,7 +85,7 @@ GridIntersection::processingRegion(const Envelope& raster_extent, const Geometry return ret; } -GridIntersection::GridIntersection(const Grid& raster_grid, const Geometry& g, std::shared_ptr> cov) +GridIntersection::GridIntersection(const Grid& raster_grid, const Geometry& g, const std::shared_ptr>& cov) : m_geometry_grid{ make_infinite(raster_grid, *g.getEnvelopeInternal()) } , m_results{ cov ? cov : std::make_shared>(raster_grid.getNumRows(), raster_grid.getNumCols()) } , m_first_geom{ true } @@ -97,9 +100,11 @@ GridIntersection::GridIntersection(const Grid& raster_grid, cons } } -GridIntersection::GridIntersection(const Grid& raster_grid, const Envelope& box, std::shared_ptr> cov) +GridIntersection::GridIntersection(const Grid& raster_grid, const Envelope& box, const std::shared_ptr>& cov) : m_geometry_grid{ make_infinite(raster_grid, box) } , m_results{ cov ? cov : std::make_shared>(raster_grid.getNumRows(), raster_grid.getNumCols()) } + , m_first_geom { true } + , m_areal{ false } { if (!m_geometry_grid.isEmpty()) { processRectangularRing(box, true); @@ -247,7 +252,7 @@ GridIntersection::setAreal(bool areal) m_areal = areal; } else { if (m_areal != areal) { - throw std::runtime_error("Mixed-type geometries not supported."); + throw util::GEOSException("Mixed-type geometries not supported."); } } } @@ -308,8 +313,8 @@ collect_lengths(const Matrix>& cells) return lengths; } -void -traverse_cells(Matrix>& cells, std::vector& coords, const Grid& ring_grid, bool areal) +static void +traverseCells(Matrix>& cells, std::vector& coords, const Grid& ring_grid, bool areal, bool exitOnBoundary, const void* parentage) { size_t pos = 0; size_t row = ring_grid.getRow(coords.front().y); @@ -323,7 +328,7 @@ traverse_cells(Matrix>& cells, std::vector& const CoordinateXY* next_coord = last_exit ? last_exit : &coords[pos]; const CoordinateXY* prev_coord = pos > 0 ? &coords[pos - 1] : nullptr; - cell.take(*next_coord, prev_coord); + cell.take(*next_coord, prev_coord, exitOnBoundary, parentage); if (cell.getLastTraversal().isExited()) { // Only push our exit coordinate if it's not same as the @@ -371,7 +376,7 @@ traverse_cells(Matrix>& cells, std::vector& col++; break; default: - throw std::runtime_error("Invalid traversal"); + throw util::GEOSException("Invalid traversal"); } } } @@ -437,7 +442,7 @@ GridIntersection::processLine(const LineString& ls, bool exterior_ring) } Matrix> cells(ring_grid.getNumRows(), ring_grid.getNumCols()); - traverse_cells(cells, coordsVec, ring_grid, m_areal); + traverseCells(cells, coordsVec, ring_grid, m_areal, false, &ls); // Compute the fraction covered for all cells and assign it to // the area matrix @@ -463,8 +468,8 @@ GridIntersection::addRingResults(size_t i0, size_t j0, const Matrix& area } } -void -traverse_ring(Matrix>& cells, const Grid& grid, const LineString& g, bool want_ccw) +static void +traverseRing(Matrix>& cells, const Grid& grid, const LineString& g, bool want_ccw, bool exitOnBoundary) { const geom::CoordinateSequence& seq = *g.getCoordinatesRO(); @@ -473,12 +478,12 @@ traverse_ring(Matrix>& cells, const Grid& if (want_ccw != algorithm::Orientation::isCCW(&seq)) { std::reverse(coords.begin(), coords.end()); - } - traverse_cells(cells, coords, grid, true); + } + traverseCells(cells, coords, grid, true, exitOnBoundary, &g); } -void -traverse_polygons(Matrix>& cells, const Grid& grid, const Geometry& g) +static void +traversePolygons(Matrix>& cells, const Grid& grid, const Geometry& g, bool exitOnBoundary) { using geom::GeometryTypeId; @@ -490,20 +495,71 @@ traverse_polygons(Matrix>& cells, const Grid(gi); const LineString& shell = *poly.getExteriorRing(); - traverse_ring(cells, grid, shell, true); + traverseRing(cells, grid, shell, true, exitOnBoundary); std::size_t nrings = poly.getNumInteriorRing(); for (std::size_t j = 0; j < nrings; j++) { - traverse_ring(cells, grid, *poly.getInteriorRingN(j), false); + traverseRing(cells, grid, *poly.getInteriorRingN(j), false, exitOnBoundary); } } else if (typ == GeometryTypeId::GEOS_MULTIPOLYGON || typ == GeometryTypeId::GEOS_GEOMETRYCOLLECTION) { - traverse_polygons(cells, grid, gi); + traversePolygons(cells, grid, gi, exitOnBoundary); } else { - throw std::runtime_error("Unsupported geometry type"); + throw util::GEOSException("Unsupported geometry type"); + } + } +} + +// Create a rectangular polygon from a set of points lying on the boundary +// of a specified envelope. The provided points do not need to include the +// corners of the envelope itself. +static std::unique_ptr +createRectangle(const geom::GeometryFactory& gfact, const Envelope& env, std::vector& points) +{ + if (points.empty()) { + return gfact.toGeometry(&env); + } else { + auto perimeterDistanceCmp = [&env](const CoordinateXY& a, const CoordinateXY& b) { + return PerimeterDistance::isLessThan(env, a, b); + }; + + points.emplace_back(env.getMinX(), env.getMinY()); + points.emplace_back(env.getMinX(), env.getMaxY()); + points.emplace_back(env.getMaxX(), env.getMinY()); + points.emplace_back(env.getMaxX(), env.getMaxY()); + + std::sort(points.begin(), points.end(), perimeterDistanceCmp); + points.push_back(points.front()); + auto seq = std::make_unique(0, false, false); + seq->reserve(points.size()); + for (const auto& c : points) { + seq->add(c, false); } + + auto ring = gfact.createLinearRing(std::move(seq)); + return gfact.createPolygon(std::move(ring)); } } +// Get any points from adjacent cells that are also on the boundary of this cell. +static std::vector +getExtraNodes(const Matrix > &cells, size_t i, size_t j) { + std::vector points; + + if (const Cell *above = cells(i - 1, j).get()) { + above->getEdgePoints(Side::BOTTOM, points); + } + if (const Cell *below = cells(i + 1, j).get()) { + below->getEdgePoints(Side::TOP, points); + } + if (const Cell *left = cells(i, j - 1).get()) { + left->getEdgePoints(Side::RIGHT, points); + } + if (const Cell *right = cells(i, j + 1).get()) { + right->getEdgePoints(Side::LEFT, points); + } + + return points; +} std::unique_ptr GridIntersection::subdividePolygon(const Grid& p_grid, const Geometry& g, bool includeExterior) @@ -512,12 +568,17 @@ GridIntersection::subdividePolygon(const Grid& p_grid, const Geo const geom::GeometryFactory& gfact = *g.getFactory(); - Grid grid = make_infinite(p_grid, *g.getEnvelopeInternal()); - Matrix> cells(grid.getNumRows(), grid.getNumCols()); + const auto cropped_grid = p_grid.shrinkToFit(p_grid.getExtent().intersection(*g.getEnvelopeInternal()), false); + + geom::Envelope gridExtent = *g.getEnvelopeInternal(); + gridExtent.expandBy(1); + + const Grid cell_grid = make_infinite(cropped_grid, gridExtent); + Matrix> cells(cell_grid.getNumRows(), cell_grid.getNumCols()); - traverse_polygons(cells, grid, g); + traversePolygons(cells, cell_grid, g, true); - const auto areas = collectAreas(cells, p_grid, g); + const auto areas = collectAreas(cells, cropped_grid, g); std::vector> geoms; std::vector> edge_geoms; @@ -529,23 +590,27 @@ GridIntersection::subdividePolygon(const Grid& p_grid, const Geo const bool col_edge = j == 0 || j == cells.getNumCols() - 1; const bool edge = row_edge || col_edge; - if (cells(i, j) != nullptr) { + const Cell* cell = cells(i, j).get(); + + if (cell != nullptr && cell->isDetermined()) { + // It is possible for the area to equal the cell area when the polygon is not the same as the + // cell area. (See GridIntersectionTest test #46). So, we only compare the area to + // fill_values::INTERIOR for cells that have not been traversed. if (edge) { if (includeExterior) { - edge_geoms.push_back(cells(i, j)->getCoveredPolygons(gfact)); + edge_geoms.push_back(cell->getCoveredPolygons(gfact)); } - } else if (areas(i - 1, j - 1) == fill_values::INTERIOR) { - // Cell is completely covered by polygon - Envelope env = cells(i, j)->box(); - geoms.push_back(gfact.toGeometry(&env)); } else if (areas(i - 1, j - 1) != fill_values::EXTERIOR) { // Cell is partly covered by polygon - geoms.push_back(cells(i, j)->getCoveredPolygons(gfact)); + geoms.push_back(cell->getCoveredPolygons(gfact)); } } else if (!edge && areas(i - 1, j - 1) == fill_values::INTERIOR) { // Cell is entirely covered by polygon - Envelope env = grid.getCellEnvelope(i, j); - geoms.push_back(gfact.toGeometry(&env)); + // In order to have the outputs forms a properly noded coverage, + // we need to add nodes from adjacent polygons. + Envelope env = cell_grid.getCellEnvelope(i, j); + std::vector points = getExtraNodes(cells, i, j); + geoms.push_back(createRectangle(gfact, env, points)); } } } @@ -553,6 +618,8 @@ GridIntersection::subdividePolygon(const Grid& p_grid, const Geo if (!edge_geoms.empty()) { auto edge_coll = gfact.createGeometryCollection(std::move(edge_geoms)); geoms.push_back(overlayng::CoverageUnion::geomunion(edge_coll.get())); + // Polygon generated here may have extra nodes that do not match adjacent polygons not processed through + // subdividePolygon. } return gfact.createGeometryCollection(std::move(geoms)); diff --git a/src/operation/grid/PerimeterDistance.cpp b/src/operation/grid/PerimeterDistance.cpp index 540680178..15c789e6f 100644 --- a/src/operation/grid/PerimeterDistance.cpp +++ b/src/operation/grid/PerimeterDistance.cpp @@ -61,4 +61,50 @@ PerimeterDistance::getPerimeterDistanceCCW(double measure1, double measure2, dou return perimeter + measure1 - measure2; } +bool +PerimeterDistance::isBetweenCCW(double x, double a, double b) +{ + if (a < b) { + return x <= a || x >= b; + } + return x <= a && x >= b; +} + +bool +PerimeterDistance::isLessThan(const geom::Envelope& env, const geom::CoordinateXY& a, const geom::CoordinateXY& b) { + const auto pda = getPerimeterDistance(env, a); + const auto pdb = getPerimeterDistance(env, b); + + if (pda < pdb) { + return true; + } + if (pda > pdb) { + return false; + } + if (a.equals2D(b)) { + return false; + } + + // Points are not equal but are so close that their perimeter distance is the same. + // Need to figure out the spatial arrangement of the points. + if (a.y == b.y) { + // top + if (a.y == env.getMaxY()) { + return a.x < b.x; + } + // bottom + if (a.y == env.getMinY()) { + return b.x > a.x; + } + } + + // left + if (a.x == env.getMinX()) { + return a.y < b.y; + } + // right + return b.y > a.y; +} + + } diff --git a/src/operation/grid/Traversal.cpp b/src/operation/grid/Traversal.cpp index d98ece5ec..035c85b53 100644 --- a/src/operation/grid/Traversal.cpp +++ b/src/operation/grid/Traversal.cpp @@ -15,6 +15,7 @@ #include #include +#include #include using geos::geom::CoordinateXY; @@ -24,7 +25,9 @@ namespace geos::operation::grid { void Traversal::add(const CoordinateXY& c) { - m_coords.push_back(c); + if (m_coords.empty() || (m_coords.back() != c)) { + m_coords.push_back(c); + } } bool @@ -34,7 +37,7 @@ Traversal::isEmpty() const } void -Traversal::enter(const CoordinateXY& c, Side s) +Traversal::enter(const CoordinateXY& c, Side s, const void* parentage) { if (!m_coords.empty()) { throw std::runtime_error("Traversal already started"); @@ -42,6 +45,7 @@ Traversal::enter(const CoordinateXY& c, Side s) add(c); m_entry = s; + m_parentage = parentage; } void @@ -57,6 +61,13 @@ Traversal::isClosedRing() const return m_coords.size() >= 3 && m_coords[0] == m_coords[m_coords.size() - 1]; } +bool +Traversal::isClosedRingWithArea() const +{ + // TODO: Don't actually need to compute the area here, just need to make sure there are >= 3 unique points. + return isClosedRing() && algorithm::Area::ofRing(m_coords) > 0; +} + bool Traversal::isEntered() const { @@ -87,10 +98,16 @@ Traversal::isTraversed() const return isEntered() && isExited(); } +const CoordinateXY& +Traversal::getFirstCoordinate() const +{ + return m_coords.front(); +} + const CoordinateXY& Traversal::getLastCoordinate() const { - return m_coords.at(m_coords.size() - 1); + return m_coords.back(); } const CoordinateXY& diff --git a/src/operation/grid/TraversalAreas.cpp b/src/operation/grid/TraversalAreas.cpp index 9b3cd59ca..9150b286a 100644 --- a/src/operation/grid/TraversalAreas.cpp +++ b/src/operation/grid/TraversalAreas.cpp @@ -13,6 +13,7 @@ **********************************************************************/ #include +#include #include #include @@ -21,10 +22,12 @@ #include #include #include +#include #include - -#include +#include #include +#include +#include using geos::geom::CoordinateSequence; using geos::geom::CoordinateXY; @@ -32,22 +35,58 @@ using geos::geom::Geometry; using geos::geom::GeometryFactory; using geos::geom::Envelope; +#define DEBUG_TRAVERSAL_AREAS 0 + +#if DEBUG_TRAVERSAL_AREAS +#include +#endif + namespace geos::operation::grid { +/// A CoordinateChain stores a set of coordinates whose start and end points lie on +/// the envelope. struct CoordinateChain { - double start; // perimeter distance value of the first coordinate - double stop; // perimeter distance value of the last coordinate - const std::vector* coordinates; + const double start; // perimeter distance value of the first coordinate + const double stop; // perimeter distance value of the last coordinate + const void* parentage; // pointer to track the linear geometry from which the coordinates originated bool visited; - CoordinateChain(double p_start, double p_stop, const std::vector* p_coords) + CoordinateChain(double p_start, double p_stop, const std::vector& p_coords) : start{ p_start } , stop{ p_stop } - , coordinates{ p_coords } + , parentage{ nullptr } , visited{ false } + , m_start{ p_coords.begin() } + , m_stop{ p_coords.end() } { } + + CoordinateChain(double p_start, double p_stop, const std::vector& p_coords, size_t from, size_t to, const void* p_parentage) + : start{ p_start } + , stop{ p_stop } + , parentage{ p_parentage } + , visited{ false } + , m_start{ std::next(p_coords.begin(), static_cast(from)) } + , m_stop{ std::next(p_coords.begin(), static_cast(to + 1)) } + { + } + + auto begin() const { + return m_start; + } + + auto end() const { + return m_stop; + } + + auto getSize() const { + return m_stop - m_start; + } + +private: + const std::vector::const_iterator m_start; + const std::vector::const_iterator m_stop; }; static double @@ -66,12 +105,30 @@ getNextChain(std::vector& chains, CoordinateChain* min = nullptr; double min_distance = std::numeric_limits::max(); +#if DEBUG_TRAVERSAL_AREAS + std::cout << "Acceptable range: " << chain->stop << " " << kill->start << std::endl; +#endif + for (CoordinateChain& candidate : chains) { if (candidate.visited && std::addressof(candidate) != kill) { continue; } - double distance = exit_to_entry_perimeter_distance_ccw(*chain, candidate, perimeter); + const double distance = exit_to_entry_perimeter_distance_ccw(*chain, candidate, perimeter); + +#if DEBUG_TRAVERSAL_AREAS + std::cout << "Distance " << distance << ", start " << candidate.start << ", stop " << candidate.stop << ": "; + for (const auto& c : candidate) { + std::cout << c << ", "; + } + std::cout<< std::endl; +#endif + + if (distance == 0 && std::addressof(candidate) != kill && !PerimeterDistance::isBetweenCCW(candidate.stop, chain->stop, kill->start)) { + // Make sure the candidate chain closes the ring instead of starting a new ring. + continue; + } + if (distance < min_distance) { min_distance = distance; min = std::addressof(candidate); @@ -98,40 +155,226 @@ hasMultipleUniqueCoordinates(const T& vec) return false; } +static bool +isRingCCW(const std::vector& coords) +{ + // TODO: Remove copy + CoordinateSequence seq(0, false, false); + seq.setPoints(coords); + return algorithm::Orientation::isCCW(&seq); +} + +static bool +isRing(const std::vector& coords) +{ + return coords.front() == coords.back() && hasMultipleUniqueCoordinates(coords); +} + +static bool +containsProperly(const Envelope& e, const std::vector& coords) +{ + return std::all_of(coords.begin(), coords.end(), [&e](const auto& c) { + return e.containsProperly(c); + }); +} + +static std::vector +scrollRingToStartOnEdge(const Envelope& box, const std::vector& coords) +{ + for (std::size_t i = 0; i < coords.size(); i++) + { + // Potentially this ring as an ordinary Traversal. + if (!box.containsProperly(coords[i])) { + // First, modify the ring so that its start point lies on the boundary. + std::vector mod; + mod.reserve(coords.size()); + + for (std::size_t j = i; j < coords.size(); j++) { + if (mod.empty() || mod.back() != coords[j]) { + mod.push_back(coords[j]); + } + } + for (std::size_t j = 0; j < i; j++) { + if (mod.empty() || mod.back() != coords[j]) { + mod.push_back(coords[j]); + } + } + mod.push_back(mod.front()); // close ring + + return mod; + } + } + + throw std::runtime_error("Ring does not touch box"); +} + +static bool +hasEdgeSegment(const Envelope& box, const std::vector& coords) +{ + for (size_t i = 1; i < coords.size(); i++) { + const CoordinateXY& p0 = coords[i - 1]; + const CoordinateXY& p1 = coords[i]; + + if (p0.x == p1.x && (p0.x == box.getMinX() || p0.x == box.getMaxX())) { + return true; + } + + if (p0.y == p1.y && (p0.y == box.getMinY() || p0.y == box.getMaxY())) { + return true; + } + } + + return false; +} + /** * @brief Identify counter-clockwise rings formed by the supplied coordinate vectors and the boundary of this box. * * @param box Box to be included in rings - * @param coord_lists A list of coordinate vectors, with each vector representing a traversal of `box` or a closed ring. + * @param traversals A list of Traversals, with each vector representing a traversal of `box` or a closed ring. * @param visitor Function be applied to each ring identified. Because `coord_lists` may include both clockwise and * counter-clockwise closed rings, `visitor` will be provided with the orientation of each ring as an * argument. */ -template +template void -visitRings(const Envelope& box, const std::vector*>& coord_lists, F&& visitor) +visitRings(const Envelope& box, const std::vector& traversals, F&& visitor) { std::vector chains; - chains.reserve(coord_lists.size() + 4); - - for (const auto& coords : coord_lists) { - if (!hasMultipleUniqueCoordinates(*coords)) { + chains.reserve(traversals.size() + 4); + std::deque> coordinateStore; + + std::vector boxPoints; // a list of nodes to be included in the boundary of the + // polygon, if it is the same as the boundary of the + // cell. Necessary to form a correctly noded polygon + // coverage. + std::optional isInterior; // whether the cell is entirely inside the polygon. + + double areaCCW = 0; // cumulative area of CCW-oriented rings + double areaCW = 0; // cumulative area of CW-oriented rings + + for (const auto& traversal : traversals) { + if (!traversal->hasMultipleUniqueCoordinates()) { + if constexpr (validPolygonsAndCoverage) { + const auto& coords = traversal->getCoordinates(); + double m = PerimeterDistance::getPerimeterDistance(box, coords.front()); + chains.emplace_back(m, m, coords); + boxPoints.push_back(coords.front()); + } continue; } - if (coords->front() == coords->back() && hasMultipleUniqueCoordinates(*coords)) { - // Closed ring. Check orientation. + const std::vector* coords = &traversal->getCoordinates(); + + if (isRing(*coords)) { + if constexpr (!validPolygonsAndCoverage) { + const bool isCCW = isRingCCW(*coords); + + visitor(*coords, isCCW, false); + + if (isCCW) { + areaCCW += algorithm::Area::ofRing(*coords); + } else { + areaCW += algorithm::Area::ofRing(*coords); + } + + continue; + } - // TODO: Remove copy - CoordinateSequence seq(0, false, false); - seq.setPoints(*coords); - bool is_ccw = algorithm::Orientation::isCCW(&seq); - visitor(*coords, is_ccw); + if (containsProperly(box, *coords) || !hasEdgeSegment(box, *coords)) { + const bool isCCW = isRingCCW(*coords); + + // Add nodes to box edges + for (const auto& c : *coords) { + if (!box.containsProperly(c)) { + boxPoints.push_back(c); + coordinateStore.emplace_back(std::vector{c}); + double m = PerimeterDistance::getPerimeterDistance(box, c); + chains.emplace_back(m, m, coordinateStore.back()); + } + } + + constexpr bool hasMultipleParents = false; + visitor(*coords, isCCW, hasMultipleParents); + + if (isCCW) { + areaCCW += algorithm::Area::ofRing(*coords); + } else { + areaCW += algorithm::Area::ofRing(*coords); + } + + continue; + } + + // TODO add node on edge. + auto modCoords = scrollRingToStartOnEdge(box, *coords); + coordinateStore.push_back(std::move(modCoords)); + coords = &coordinateStore.back(); + } + + if constexpr (validPolygonsAndCoverage) { + // Split coordinates into separate chains when they touch an edge + // This prevents the creation of self-touching rings. + // For area calculations, this doesn't matter, and the step can be skipped. + size_t from = 0; + size_t uniquePoints = 1; + + for (size_t to = 1; to < coords->size(); to++) { // cppcheck-suppress invalidContainer + const CoordinateXY& prev = (*coords)[from]; + const CoordinateXY& c = (*coords)[to]; + + if (!c.equals2D(prev)) { + uniquePoints++; + } + + const bool ptOnTop = c.y == box.getMaxY(); + const bool ptOnBottom = c.y == box.getMinY(); + const bool ptOnLeft = c.x == box.getMinX(); + const bool ptOnRight = c.x == box.getMaxX(); + + if (ptOnTop || ptOnBottom || ptOnLeft || ptOnRight) { + bool isEdgeSegment = false; + if (uniquePoints == 2) { + if (ptOnTop && prev.y == c.y) { + isEdgeSegment = true; + isInterior = c.x < prev.x; + } else if (ptOnBottom && prev.y == c.y) { + isEdgeSegment = true; + isInterior = c.x > prev.x; + } else if (ptOnLeft && prev.x == c.x) { + isEdgeSegment = true; + isInterior = c.y < prev.y; + } else if (ptOnRight && prev.x == c.x) { + isEdgeSegment = true; + isInterior = c.y > prev.y; + } + } + + if (isEdgeSegment) { + boxPoints.push_back(c); + boxPoints.push_back(prev); + + coordinateStore.emplace_back(std::vector{c}); + double m = PerimeterDistance::getPerimeterDistance(box, c); + chains.emplace_back(m, m, coordinateStore.back()); + + coordinateStore.emplace_back(std::vector{prev}); + m = PerimeterDistance::getPerimeterDistance(box, prev); + chains.emplace_back(m, m, coordinateStore.back()); + + } else { + double start = PerimeterDistance::getPerimeterDistance(box, (*coords)[from]); + double stop = PerimeterDistance::getPerimeterDistance(box, (*coords)[to]); + chains.emplace_back(start, stop, *coords, from, to, traversal->getParentage()); + } + from = to; + uniquePoints = 1; + } + } } else { double start = PerimeterDistance::getPerimeterDistance(box, coords->front()); double stop = PerimeterDistance::getPerimeterDistance(box, coords->back()); - - chains.emplace_back(start, stop, coords); + chains.emplace_back(start, stop, *coords); } } @@ -146,113 +389,186 @@ visitRings(const Envelope& box, const std::vector bottom_right = { CoordinateXY(box.getMaxX(), box.getMinY()) }; // Add chains for corners - chains.emplace_back(0.0, 0.0, &bottom_left); - chains.emplace_back(height, height, &top_left); - chains.emplace_back(height + width, height + width, &top_right); - chains.emplace_back(2 * height + width, 2 * height + width, &bottom_right); + chains.emplace_back(0.0, 0.0, bottom_left); + chains.emplace_back(height, height, top_left); + chains.emplace_back(height + width, height + width, top_right); + chains.emplace_back(2 * height + width, 2 * height + width, bottom_right); std::vector coords; + +#if DEBUG_TRAVERSAL_AREAS + std::cout << std::endl << std::fixed; + std::cout << std::setprecision(18); + std::cout << "Identifying rings in box " << box.getMinX() << ":" << box.getMaxX() << "," << box.getMinY() << ":" << box.getMaxY() << std::endl; + std::cout << std::setprecision(18); + std::cout << "Available chains:" << std::endl; + for (const auto& chain : chains) { + for (const auto& coord : chain) { + std:: cout << coord << ", "; + } + std::cout << std::endl; + } +#endif + for (auto& chain_ref : chains) { - if (chain_ref.visited || chain_ref.coordinates->size() == 1) { + if (chain_ref.visited || chain_ref.getSize() == 1) { continue; } coords.clear(); - +#if DEBUG_TRAVERSAL_AREAS + std::cout << "New ring." << std::endl; +#endif CoordinateChain* chain = std::addressof(chain_ref); - CoordinateChain* first_chain = chain; + const CoordinateChain* first_chain = chain; + bool hasMultipleParents = false; + constexpr bool isCCW = true; do { chain->visited = true; - coords.insert(coords.end(), chain->coordinates->cbegin(), chain->coordinates->cend()); + coords.insert(coords.end(), chain->begin(), chain->end()); + + if (chain->parentage != nullptr && chain->parentage != first_chain->parentage) { + hasMultipleParents = true; + } + +#if DEBUG_TRAVERSAL_AREAS + std::cout << "Added chain "; + for (const auto& c : *chain) { + std::cout << c << ", "; + } + std::cout<< std::endl; +#endif + chain = getNextChain(chains, chain, first_chain, perimeter); } while (chain != first_chain); coords.push_back(coords[0]); +#if DEBUG_TRAVERSAL_AREAS + std::cout << "Completed chain LINESTRING ("; + { + bool first = true; + for (const auto& c : coords) { + if (first) { + first = false; + } else { + std::cout << ", "; + } + std::cout << c; + } + std::cout << ")"; + } + std::cout<< std::endl; +#endif if (hasMultipleUniqueCoordinates(coords)) { - visitor(coords, true); + visitor(coords, isCCW, hasMultipleParents); + double area = algorithm::Area::ofRing(coords); + areaCCW += area; + if (area == 0) { + isInterior = false; + } } } + + if (areaCCW == 0 && areaCW == 0 && !isInterior.has_value()) { + throw util::GEOSException("Cannot determine coverage fraction (it is either 0 or 100%)"); + } + + if (areaCW > areaCCW || (areaCCW == 0 && areaCW == 0 && isInterior.value())) { + // Ring is the same as the cell envelope + const bool sortIsNeeded = !boxPoints.empty(); + + boxPoints.emplace_back(box.getMinX(), box.getMinY()); + boxPoints.emplace_back(box.getMaxX(), box.getMinY()); + boxPoints.emplace_back(box.getMaxX(), box.getMaxY()); + boxPoints.emplace_back(box.getMinX(), box.getMaxY()); + + if (sortIsNeeded) { + std::sort(boxPoints.begin(), boxPoints.end(), [&box](const CoordinateXY& lhs, const CoordinateXY& rhs) { + return PerimeterDistance::isLessThan(box, lhs, rhs); + }); + } + + boxPoints.push_back(boxPoints.front()); + + visitor(boxPoints, true, false); + } } double -TraversalAreas::getLeftHandArea(const Envelope& box, const std::vector*>& coord_lists) +TraversalAreas::getLeftHandArea(const Envelope& box, const std::vector& coord_lists) { double ccw_sum = 0; double cw_sum = 0; - bool found_a_ring = false; - visitRings(box, coord_lists, [&cw_sum, &ccw_sum, &found_a_ring](const std::vector& coords, bool is_ccw) { - found_a_ring = true; + constexpr bool validPolygons = false; - if (is_ccw) { + visitRings(box, coord_lists, [&cw_sum, &ccw_sum](const std::vector& coords, bool isCCW, bool) { + if (isCCW) { ccw_sum += algorithm::Area::ofRing(coords); } else { cw_sum += algorithm::Area::ofRing(coords); } }); - if (!found_a_ring) { - throw std::runtime_error("Cannot determine coverage fraction (it is either 0 or 100%)"); - } - - // If this box has only clockwise rings (holes) then the area - // of those holes should be subtracted from the complete box area. - if (ccw_sum < cw_sum) { - ccw_sum += box.getArea(); - } - return ccw_sum - cw_sum; } std::unique_ptr -TraversalAreas::getLeftHandRings(const GeometryFactory& gfact, const Envelope& box, const std::vector*>& coord_lists) -{ +TraversalAreas::getLeftHandRings(const GeometryFactory& gfact, const Envelope& box, const std::vector& coord_lists) { using geom::LinearRing; std::vector> shells; std::vector> holes; - bool found_a_ring = false; + bool checkValidity = false; + + constexpr bool validPolygons = true; - visitRings(box, coord_lists, [&gfact, &shells, &holes, &found_a_ring](const std::vector& coords, bool is_ccw) { - found_a_ring = true; + visitRings(box, coord_lists, [&gfact, &shells, &holes, &checkValidity](const std::vector& coords, bool isCCW, bool hasMultipleParents) { + // If a given ring was created from traversals from both a ring and shell, for example, it is possible for + // the ring to self-intersect. Rather than try and detect self-intersections on the fly (rare?) we check and + // repair validity in this limited case. + checkValidity |= hasMultipleParents; - // finding a collapsed ring is sufficient to determine whether the cell interior is inside our outside, + // finding a collapsed ring is sufficient to determine whether the cell interior is inside or outside, // but we don't want to actually construct the ring. if (algorithm::Area::ofRing(coords) == 0) { return; } - CoordinateSequence seq(0, false, false); - seq.setPoints(coords); + auto seq = std::make_unique(0, false, false); + seq->reserve(coords.size()); + seq->setPoints(coords); + if (seq->hasRepeatedPoints()) { + seq = valid::RepeatedPointRemover::removeRepeatedPoints(seq.get()); + } auto ring = gfact.createLinearRing(std::move(seq)); - if (is_ccw) { + if (isCCW) { shells.push_back(std::move(ring)); } else { holes.push_back(std::move(ring)); } + }); - if (!found_a_ring) { - throw std::runtime_error("Cannot determine coverage fraction (it is either 0 or 100%)"); +#if DEBUG_TRAVERSAL_AREAS + std::cout << "SHELLS" << std::endl; + for (const auto& shell: shells) { + std::cout << shell->toString() << std::endl; + } + std::cout << "HOLES" << std::endl; + for (const auto& hole: holes) { + std::cout << hole->toString() << std::endl; } +#endif if (shells.empty() && holes.empty()) { return gfact.createPolygon(); } - if (shells.empty() && !holes.empty()) { - CoordinateSequence seq(5, false, false); - seq.setAt(CoordinateXY{box.getMinX(), box.getMinY()}, 0); - seq.setAt(CoordinateXY{box.getMaxX(), box.getMinY()}, 1); - seq.setAt(CoordinateXY{box.getMaxX(), box.getMaxY()}, 2); - seq.setAt(CoordinateXY{box.getMinX(), box.getMaxY()}, 3); - seq.setAt(CoordinateXY{box.getMinX(), box.getMinY()}, 4); - - shells.push_back(gfact.createLinearRing(std::move(seq))); - } + std::unique_ptr result; if (holes.empty()) { std::vector> polygons; @@ -261,27 +577,44 @@ TraversalAreas::getLeftHandRings(const GeometryFactory& gfact, const Envelope& b } if (polygons.size() == 1) { - return std::move(polygons.front()); + result = std::move(polygons.front()); } else { - return gfact.createMultiPolygon(std::move(polygons)); + result = gfact.createMultiPolygon(std::move(polygons)); + } + } else if (shells.size() == 1) { + result = gfact.createPolygon(std::move(shells.front()), std::move(holes)); + } else { + polygonize::Polygonizer polygonizer(true); + + std::vector> holder; + + for (const auto& shell : shells) { + if (checkValidity && !shell->isSimple()) { + holder.push_back(shell->Union()); + polygonizer.add(holder.back().get()); + } else { + polygonizer.add(static_cast(shell.get())); + } + } + for (const auto& hole : holes) { + if (checkValidity && !hole->isSimple()) { + holder.push_back(hole->Union()); + polygonizer.add(holder.back().get()); + } else { + polygonizer.add(static_cast(hole.get())); + } } - } - if (shells.size() == 1) { - return gfact.createPolygon(std::move(shells.front()), std::move(holes)); + auto polygonized = polygonizer.getPolygons(); + result = gfact.createMultiPolygon(std::move(polygonized)); + checkValidity = false; } - polygonize::Polygonizer polygonizer(true); - - for (const auto& shell : shells) { - polygonizer.add(static_cast(shell.get())); - } - for (const auto& hole : holes) { - polygonizer.add(static_cast(hole.get())); + if (checkValidity && !result->isValid()) { + result = geom::util::GeometryFixer::fix(result.get()); } - auto polygonized = polygonizer.getPolygons(); - return gfact.createMultiPolygon(std::move(polygonized)); + return result; } } diff --git a/src/operation/union/OverlapUnion.cpp b/src/operation/union/OverlapUnion.cpp index dbc638772..aef18801b 100644 --- a/src/operation/union/OverlapUnion.cpp +++ b/src/operation/union/OverlapUnion.cpp @@ -196,22 +196,11 @@ intersects(const Envelope& env, const Coordinate& p0, const Coordinate& p1) return env.intersects(p0) || env.intersects(p1); } -/* static */ -static bool -containsProperly(const Envelope& env, const Coordinate& p) -{ - if (env.isNull()) return false; - return p.x > env.getMinX() && - p.x < env.getMaxX() && - p.y > env.getMinY() && - p.y < env.getMaxY(); -} - /* static */ static bool containsProperly(const Envelope& env, const Coordinate& p0, const Coordinate& p1) { - return containsProperly(env, p0) && containsProperly(env, p1); + return env.containsProperly(p0) && env.containsProperly(p1); } /* private */ diff --git a/tests/unit/capi/GEOSSubdivideByGridTest.cpp b/tests/unit/capi/GEOSSubdivideByGridTest.cpp new file mode 100644 index 000000000..d734c5bca --- /dev/null +++ b/tests/unit/capi/GEOSSubdivideByGridTest.cpp @@ -0,0 +1,56 @@ +#include +#include + +#include "capi_test_utils.h" + +namespace tut { + +struct test_capigeossubdividebygrid_data : public capitest::utility {}; + +typedef test_group group; +typedef group::object object; + +group test_capigeossubdividebygrid_group("capi::GEOSSubdivideByGrid"); + +template<> +template<> +void object::test<1>() +{ + set_test_name("rectangle overlapping grid"); + + input_ = fromWKT("POLYGON ((0.5 0.5, 2.5 0.5, 2.5 2.5, 0.5 2.5, 0.5 0.5))"); + expected_ = fromWKT("GEOMETRYCOLLECTION (" + "POLYGON ((2 2.5, 1 2.5, 1 2, 2 2, 2 2.5))," + "POLYGON ((2.5 2, 2.5 2.5, 2 2.5, 2 2, 2.5 2))," + "POLYGON ((1 1, 2 1, 2 2, 1 2, 1 1))," + "POLYGON ((2.5 1, 2.5 2, 2 2, 2 1, 2.5 1))," + "POLYGON ((1 0.5, 2 0.5, 2 1, 1 1, 1 0.5))," + "POLYGON ((2 0.5, 2.5 0.5, 2.5 1, 2 1, 2 0.5)))"); + + result_ = GEOSSubdivideByGrid(input_, 1, 0, 5, 3, 4, 3, false); + + ensure_geometry_equals_identical(expected_, result_); +} + +template<> +template<> +void object::test<2>() +{ + set_test_name("includeExterior = true"); + + input_ = fromWKT("POLYGON ((0.5 0.5, 2.5 0.5, 2.5 2.5, 0.5 2.5, 0.5 0.5))"); + expected_ = fromWKT("GEOMETRYCOLLECTION (" + "POLYGON ((2 2.5, 1 2.5, 1 2, 2 2, 2 2.5))," + "POLYGON ((2.5 2, 2.5 2.5, 2 2.5, 2 2, 2.5 2))," + "POLYGON ((1 1, 2 1, 2 2, 1 2, 1 1))," + "POLYGON ((2.5 1, 2.5 2, 2 2, 2 1, 2.5 1))," + "POLYGON ((1 0.5, 2 0.5, 2 1, 1 1, 1 0.5))," + "POLYGON ((2 0.5, 2.5 0.5, 2.5 1, 2 1, 2 0.5))," + "POLYGON ((0.5 2.5, 1 2.5, 1 2, 1 1, 1 0.5, 0.5 0.5, 0.5 1, 0.5 2, 0.5 2.5)))"); + + result_ = GEOSSubdivideByGrid(input_, 1, 0, 5, 3, 4, 3, true); + + ensure_geometry_equals_identical(expected_, result_); +} + +} \ No newline at end of file diff --git a/tests/unit/geom/EnvelopeTest.cpp b/tests/unit/geom/EnvelopeTest.cpp index b1250a63e..a816f84f3 100644 --- a/tests/unit/geom/EnvelopeTest.cpp +++ b/tests/unit/geom/EnvelopeTest.cpp @@ -655,5 +655,23 @@ void object::test<20> ensure_equals(set.size(), 2u); } -} // namespace tut +template<> +template<> +void object::test<21>() +{ + set_test_name("containsProperly"); + + Envelope env; + ensure(env.isNull()); + + ensure(!env.containsProperly({0, 0})); + env.expandToInclude({3, 7}); + env.expandToInclude( {5, 9}); + + ensure(!env.containsProperly({3, 7})); + ensure(!env.containsProperly({3, 8})); + ensure(env.containsProperly({3.0000001, 8})); +} + +} diff --git a/tests/unit/operation/grid/GridIntersectionTest.cpp b/tests/unit/operation/grid/GridIntersectionTest.cpp index c7bf4965b..2c5cb8a01 100644 --- a/tests/unit/operation/grid/GridIntersectionTest.cpp +++ b/tests/unit/operation/grid/GridIntersectionTest.cpp @@ -1,11 +1,14 @@ #include +#include #include +#include #include #include #include + using namespace geos::operation::grid; using geos::geom::CoordinateXY; using geos::geom::Envelope; @@ -14,6 +17,7 @@ using geos::geom::Geometry; namespace tut { struct test_gridintersectiontest_data : GEOSTestBase { geos::io::WKTReader wkt_reader_; + geos::io::WKBReader wkb_reader_; static void check_cell_intersections(const Matrix& actual, const std::vector>& v) @@ -22,6 +26,47 @@ struct test_gridintersectiontest_data : GEOSTestBase { ensure_equals(actual, expected); } + + static void + check_area(const Matrix& actual, const Grid ext, const Geometry& input) + { + double tot_area = 0; + for (float f : actual) { + tot_area += static_cast(f) * ext.dx()*ext.dy(); + } + + ensure_equals("Area does not equal input", tot_area, input.getArea(), 1e-6*tot_area); + } + + static void + check_subdivided_polygon(const Geometry& input, const Geometry& subdivided) + { + double tot_area = 0; + + std::vector components; + + for (size_t i = 0; i < subdivided.getNumGeometries(); i++) { + const Geometry* subg = subdivided.getGeometryN(i); + components.push_back(subg); + + ensure("subdivided component " + subg->toString() + " is invalid",subg->isValid()); + tot_area += subg->getArea(); + } + + if (!geos::coverage::CoverageValidator::isValid(components)) + { + auto invalidEdges = geos::coverage::CoverageValidator::validate(components); + invalidEdges.erase(std::remove_if(invalidEdges.begin(), invalidEdges.end(), [](const auto& edge) { + return edge == nullptr; + }), invalidEdges.end()); + auto invalidEdgeGeom = input.getFactory()->buildGeometry(std::move(invalidEdges)); + std::string message = "subdivided polygons do not form a valid coverage.\nsubdivided: " + subdivided.toString() + "\ninvalid edges: " + invalidEdgeGeom->toString(); + fail(message); + } + + std::string error = "subdivided polygon area does not match input: " + subdivided.toString(); + ensure_equals(error, tot_area, input.getArea(), input.getArea() * 1e-6); + } }; typedef test_group group; @@ -445,8 +490,11 @@ void object::test<23>() Grid extent{ { -180.5, 180.5, -90.5, 90.5 }, 0.5, 0.5 }; auto g = wkt_reader_.read("MULTIPOLYGON (((178.3736000000001 -17.33992000000002, 178.71806000000007 -17.62845999999996, 178.5527099999999 -18.150590000000008, 177.93266000000008 -18.287990000000036, 177.38145999999992 -18.164319999999975, 177.28504000000007 -17.72464999999997, 177.67087 -17.381139999999974, 178.12557000000007 -17.50480999999995, 178.3736000000001 -17.33992000000002)), ((179.36414266196417 -16.801354076946836, 178.7250593629972 -17.012041674368007, 178.5968385951172 -16.63915000000003, 179.0966093629972 -16.43398427754741, 179.4135093629972 -16.379054277547382, 180.00000000000003 -16.06713266364241, 180.00000000000003 -16.555216566639146, 179.36414266196417 -16.801354076946836)), ((-179.91736938476527 -16.501783135649347, -179.99999999999997 -16.555216566639146, -179.99999999999997 -16.06713266364241, -179.79332010904858 -16.020882256741217, -179.91736938476527 -16.501783135649347)))"); + auto rci = GridIntersection::getIntersectionFractions(extent, *g); + check_area(*rci, extent, *g); - ensure_NO_THROW(GridIntersection(extent, *g)); + auto subd = GridIntersection::subdividePolygon(extent, *g, false); + check_subdivided_polygon(*g, *subd); } template<> @@ -490,14 +538,20 @@ template<> void object::test<26>() { set_test_name("Robustness regression test #1"); + skip_on_ci(); // This test exercises some challenging behavior where a polygon follows // ymin, but the grid resolution is such that ymin < (ymax - ny*dy) - Grid extent{ { -180, 180, -90, 90 }, 1.0 / 6, 1.0 / 6 }; + Grid ext{ { -180, 180, -90, 90 }, 1.0 / 6, 1.0 / 6 }; auto g = wkt_reader_.read(load_resource("antarctica.wkt")); - ensure_NO_THROW(GridIntersection(extent, *g)); + + auto rci = GridIntersection::getIntersectionFractions(ext, *g); + check_area(*rci, ext, *g); + + auto subd = GridIntersection::subdividePolygon(ext, *g, false); + check_subdivided_polygon(*g, *subd); } template<> @@ -505,15 +559,19 @@ template<> void object::test<27>() { set_test_name("Robustness regression test #2"); + skip_on_ci(); // This test exercises some challenging behavior where a polygon follows // xmax, but the grid resolution is such that xmax < (xmin + nx*m_dx) - - Grid extent{ { -180, 180, -90, 90 }, 1.0 / 6, 1.0 / 6 }; + Grid ext{ { -180, 180, -90, 90 }, 1.0 / 6, 1.0 / 6 }; auto g = wkt_reader_.read(load_resource("russia.wkt")); - ensure_NO_THROW(GridIntersection(extent, *g)); + auto rci = GridIntersection::getIntersectionFractions(ext, *g); + check_area(*rci, ext, *g); + + auto subd = GridIntersection::subdividePolygon(ext, *g, false); + check_subdivided_polygon(*g, *subd); } template<> @@ -529,7 +587,11 @@ void object::test<29>() Envelope env = g->getEnvelopeInternal()->intersection(extent.getExtent()); extent = extent.shrinkToFit(env); - ensure_NO_THROW( GridIntersection::getIntersectionFractions(extent, *g)); + auto rci = GridIntersection::getIntersectionFractions(extent, *g); + check_area(*rci, extent, *g); + + auto subd = GridIntersection::subdividePolygon(extent, *g, false); + check_subdivided_polygon(*g, *subd); } template<> @@ -537,12 +599,17 @@ template<> void object::test<30>() { set_test_name("Robustness regression test #4"); + skip_on_ci(); - Grid extent{ { -166.84166666666667, -152.625, 66.991666666666674, 71.358333333333334 }, 0.0083333333333333332, 0.0083333333333333332 }; + Grid ext{ { -166.84166666666667, -152.625, 66.991666666666674, 71.358333333333334 }, 0.0083333333333333332, 0.0083333333333333332 }; auto g = wkt_reader_.read(load_resource("regression4.wkt")); - ensure_NO_THROW(GridIntersection(extent, *g)); + auto rci = GridIntersection::getIntersectionFractions(ext, *g); + check_area(*rci, ext, *g); + + auto subd = GridIntersection::subdividePolygon(ext, *g, false); + check_subdivided_polygon(*g, *subd); } template<> @@ -551,11 +618,15 @@ void object::test<31>() { set_test_name("robustness regression test #5"); - Grid extent{ { 0, 10, 0, 10 }, 1, 1 }; + Grid ext{ { 0, 10, 0, 10 }, 1, 1 }; auto g = wkt_reader_.read("POINT (2 2)")->buffer(1, 30); - ensure_NO_THROW( GridIntersection::getIntersectionFractions(extent, *g)); + auto rci = GridIntersection::getIntersectionFractions(ext, *g); + check_area(*rci, ext, *g); + + auto subd = GridIntersection::subdividePolygon(ext, *g, false); + check_subdivided_polygon(*g, *subd); } template<> @@ -564,24 +635,15 @@ void object::test<32>() { set_test_name("Robustness regression test #6"); - Grid ex{ { 145.925, 147.375, -35.525, -33.475 }, 0.05, 0.05 }; + Grid ext{ { 145.925, 147.375, -35.525, -33.475 }, 0.05, 0.05 }; auto g = wkt_reader_.read(load_resource("regression6.wkt")); - GridIntersection gi(ex, *g); - const auto& result = *gi.getResults(); - float tot = 0; + auto rci = GridIntersection::getIntersectionFractions(ext, *g); + check_area(*rci, ext, *g); - for (size_t i = 0; i < result.getNumRows(); i++) { - for (size_t j = 0; j < result.getNumCols(); j++) { - tot += result(i, j); - if (result(i, j) < 0 || result(i, j) > 1) { - fail(); - } - } - } - - ensure_equals(tot, 823.0f); + auto subd = GridIntersection::subdividePolygon(ext, *g, false); + check_subdivided_polygon(*g, *subd); } template<> @@ -592,6 +654,7 @@ void object::test<33>() Grid ex{ { 487800, 492800, 5813800, 5818800 }, 100, 100 }; + // Polygon is a tiny sliver with an area of 1e-8 auto g = wkt_reader_.read("POLYGON ((492094.9283999996 5816959.8553, 492374.9335527361 5816811.352641133, 492374.9335527363 5816811.352641133, 492094.9283999996 5816959.8553))"); ex = ex.shrinkToFit(*g->getEnvelopeInternal()); @@ -611,7 +674,8 @@ void object::test<33>() template<> template<> -void object::test<34>() {set_test_name("Processing region is empty when there are no polygons"); +void object::test<34>() { + set_test_name("Processing region is empty when there are no polygons"); Envelope raster_extent{ 0, 10, 0, 10 }; @@ -645,8 +709,7 @@ void object::test<36>() auto g = wkt_reader_.read("POLYGON ((1 3, 9 5, 8 9, 1 3))") ; auto subdivided = GridIntersection::subdividePolygon(ext, *g, false); - - ensure_equals(g->getArea(), subdivided->getArea()); + check_subdivided_polygon(*g, *subdivided); } template<> @@ -661,8 +724,7 @@ void object::test<37>() auto g = wkt_reader_.read("POLYGON ((8.5 8.7, 12 8, 12 12, 8 12, 8.5 8.7))"); auto subdivided = GridIntersection::subdividePolygon(ext, *g, true); - - ensure_equals("", subdivided->getArea(), g->getArea(), 1e-8); + check_subdivided_polygon(*g, *subdivided); } template<> @@ -676,9 +738,8 @@ void object::test<38>() auto g = wkt_reader_.read("POINT (5 5)")->buffer(20); - auto subdivided = GridIntersection::subdividePolygon(ext, *g, true); - - ensure_equals("", subdivided->getArea(), g->getArea(), 1e-8); + auto subd = GridIntersection::subdividePolygon(ext, *g, true); + check_subdivided_polygon(*g, *subd); } template<> @@ -692,9 +753,8 @@ void object::test<39>() auto g = GeometryFactory::getDefaultInstance()->toGeometry(&e); - auto subdivided = GridIntersection::subdividePolygon(ext, *g, false); - - ensure_equals("", subdivided->getArea(), e.getArea(), 1e-8); + auto subd = GridIntersection::subdividePolygon(ext, *g, false); + check_subdivided_polygon(*g, *subd); } template<> @@ -708,9 +768,8 @@ void object::test<40>() auto g = wkt_reader_.read("MULTIPOLYGON (((1 1, 15 1, 15 25, 1 25, 1 1), (12 12, 12 14, 14 14, 14 12, 12 12)), ((16 1, 25 1, 25 25, 16 25, 16 1)))"); - auto subdivided = GridIntersection::subdividePolygon(ext, *g, false); - - ensure_equals("", subdivided->getArea(), g->getArea(), 1e-8); + auto subd = GridIntersection::subdividePolygon(ext, *g, false); + check_subdivided_polygon(*g, *subd); } template<> @@ -724,14 +783,11 @@ void object::test<41>() auto g = wkt_reader_.read("POLYGON ((5 5, 25 5, 25 25, 5 25, 5 5), (12 12, 12 14, 14 14, 12 12))"); - auto fracs = GridIntersection::getIntersectionFractions(ext, *g); - - float sum = 0; - for (const auto& frac : *fracs) { - sum += frac * 100; - } + auto rci = GridIntersection::getIntersectionFractions(ext, *g); + check_area(*rci, ext, *g); - ensure_equals("", static_cast(sum), g->getArea(), 1e-5); + auto subd = GridIntersection::subdividePolygon(ext, *g, false); + check_subdivided_polygon(*g, *subd); } template<> @@ -779,4 +835,167 @@ void object::test<44>() ensure_THROW(GridIntersection::getIntersectionFractions(ext, *g), std::exception); } +template<> +template<> +void object::test<45>() +{ + set_test_name("subdivide polygon whose edges follow cell boundaries"); + + Envelope e(0, 10, 0, 10); + Grid ext(e, 1, 1); + + auto g = wkt_reader_.read("POLYGON ((4 0, 6 0, 6 2, 8 2, 6 4, 8 4, 5 7, 2 4, 4 4, 2 2, 4 2, 4 0))"); + auto rci = GridIntersection::getIntersectionFractions(ext, *g); + check_area(*rci, ext, *g); + + auto subd = GridIntersection::subdividePolygon(ext, *g, false); + check_subdivided_polygon(*g, *subd); +} + +template<> +template<> +void object::test<46>() +{ + set_test_name("subdivide polygon whose edges follow cell boundaries (2)"); + + Envelope e(0, 10, 0, 10); + Grid ext(e, 1, 1); + + auto g = wkt_reader_.read("POLYGON ((4.5 0, 6.5 0, 6.5 2, 8.5 2, 6.5 4, 8.5 4, 5.5 7, 2.5 4, 4.5 4, 2.5 2, 4.5 2, 4.5 0))"); + + auto rci = GridIntersection::getIntersectionFractions(ext, *g); + check_area(*rci, ext, *g); + + auto subd = GridIntersection::subdividePolygon(ext, *g, false); + check_subdivided_polygon(*g, *subd); +} + +template<> +template<> +void object::test<47>() +{ + set_test_name("subdivide polygon whose edges follow cell boundaries (3)"); + + Envelope e(0, 10, 0, 10); + Grid ext(e, 1, 1); + + auto g = wkt_reader_.read("POLYGON ((4.5 0, 6.5 0, 6.5 2, 6.8 2, 6.5 4, 6.8 4, 5.5 7, 4.2 4, 4.5 4, 4.2 2, 4.5 2, 4.5 0))"); + + auto rci = GridIntersection::getIntersectionFractions(ext, *g); + check_area(*rci, ext, *g); + + auto subd = GridIntersection::subdividePolygon(ext, *g, false); + check_subdivided_polygon(*g, *subd); +} + +template<> +template<> +void object::test<48>() +{ + set_test_name("valid polygon coverage obtained when a traversed cell covered area ~= cell area"); + + Envelope e(-180, 180, -90, 90); + Grid ext(e, 0.5, 0.5); + + auto g = wkt_reader_.read("Polygon ((-179.99999999999991473 70.99201035500004764, -179.99999999999991473 71.05263157900003534, -179.99999999999991473 71.5366879880000397, -179.86286373599992316 71.53864166900007149, -179.91222083199991744 71.55585358300004373, -179.90074622299991347 71.55849844000005078, -179.79881751199991413 71.56907786700003271, -179.75743567599991479 71.58319733300004373, -179.73595130099991479 71.58641185100003668, -179.7154434889999095 71.58323802300003535, -179.69749915299991017 71.57733795800004373, -179.67870032499990884 71.57367584800005034, -179.61082923099991149 71.58519114800003535, -179.37205969999990884 71.56907786700003271, -179.3267716139999095 71.55548737200007281, -179.30683346299991854 71.55756256700004769, -179.28718014199992581 71.56293366100004505, -179.24286861899992118 71.56907786700003271, -179.20466061099992316 71.58319733300004373, -179.07457434799991347 71.60004303600004505, -178.73471025299991766 71.57037995000004571, -178.39484615799992184 71.54071686400004637, -178.32319088399989937 71.51837799700007281, -178.25963294199991083 71.51068756700004769, -178.30488033799991854 71.51312897300005034, -178.32347571499991545 71.51512278900003139, -178.3415421209999181 71.51752350500004241, -178.32245846299991854 71.50543854400007149, -178.21532141799991678 71.47801341400003139, -178.19347083199991744 71.47662995000007413, -178.14777584499989871 71.48517487200007281, -178.12446041599991986 71.48187897300005034, -178.00572669199991083 71.44863515800005871, -178.01720130099991479 71.44139232000003403, -178.05418860599991149 71.42877838700007942, -178.04706783799991854 71.42572663000004241, -178.03343665299991017 71.4177920590000781, -178.02623450399991611 71.41510651200007942, -178.03010006399992449 71.41347890800005871, -178.03990637899991611 71.40766022300005034, -177.97089596299991854 71.39642975500004241, -177.77985592399991788 71.33319733300004373, -177.71837317599991479 71.30524323100007678, -177.70641028599990818 71.30390045800004373, -177.68211829299991678 71.30487702000004901, -177.67027747299991347 71.30182526200007942, -177.65538489499991215 71.29315827000004901, -177.58759518099992647 71.28595612200007281, -177.5485733709999181 71.29486725500004241, -177.53111731699991083 71.29633209800005034, -177.51410885299992515 71.29340241100004505, -177.4986466139999095 71.28473541900007149, -177.50621497299991347 71.26862213700007942, -177.48700924399992118 71.25873444200004769, -177.45970618399991281 71.24990469000005078, -177.44343014199992581 71.23700592700004108, -177.4459122389999095 71.22264232000003403, -177.45775305899991281 71.20937734600005342, -177.50780188699991413 71.17377350500004241, -177.58116614499991215 71.1476097680000521, -177.63764400899989937 71.1170108090000781, -177.68415279899991788 71.11098867400005474, -177.7519018219999225 71.09296295800004373, -177.81928463399989937 71.08466217700004108, -177.87767493399991281 71.0525576840000781, -177.93049068899992449 71.04144928600004505, -178.20661373599992316 71.03839752800007545, -178.31012936099992316 71.01361725500004241, -178.59302730999991127 70.99732086800005959, -178.87592525899989937 70.98102448100007678, -178.9802953769999192 70.95066966400003139, -179.34211178299992184 70.9080264340000781, -179.33625240799992184 70.91107819200004769, -179.32225501199991413 70.9216983090000781, -179.36449133999991545 70.93024323100007678, -179.45750891799991678 70.91551341400003139, -179.50121008999991545 70.919663804000038, -179.66600501199991413 70.96548086100006003, -179.85338294199991083 70.97943756700004769, -179.88878333199991744 70.99359772300005034, -179.90754146999989871 70.99677155200004108, -179.99999999999991473 70.99201035500004764))"); + + auto subd = GridIntersection::subdividePolygon(ext, *g, false); + + check_subdivided_polygon(*g, *subd); +} + +template<> +template<> +void object::test<49>() +{ + set_test_name("self-touching rings force geometry to be corrected"); + + Envelope e(3000000, 10000000, 525000, 6595000); + Grid ext(e, 10000, 10000); + + auto g = wkt_reader_.read("MultiPolygon (((5196000 2052000, 5184185 2054473, 5182537 2054890, 5182796 2055916, 5182006 2056057, 5182183 2056774, 5181023 2058767, 5180374 2058127, 5180034 2058226, 5179989 2057895, 5179854 2057364, 5179674 2056658, 5179236 2056764, 5179289 2055146, 5180169 2052000, 5175958 2052000, 5175900 2068000, 5196000 2068000, 5196000 2052000),(5183832 2056356, 5183529 2055571, 5184300 2055372, 5184506 2056186, 5183832 2056356),(5179491 2062463, 5179636 2059879, 5180441 2059638, 5180438 2059666, 5180855 2061320, 5180076 2061534, 5180260 2062270, 5179491 2062463),(5181043 2062069, 5181035 2065286, 5180449.50561340618878603 2064531.85610372573137283, 5179685 2063225, 5180476 2063057, 5180260 2062270, 5181043 2062069)),((5180501 2056431, 5180341 2055640, 5179465 2055835, 5179674 2056658, 5180501 2056431)),((5180501 2056431, 5180641 2057164, 5181406 2056969, 5181239 2056302, 5180501 2056431))) "); + auto subd = GridIntersection::subdividePolygon(ext, *g, false); + + check_subdivided_polygon(*g, *subd); +} + +template<> +template<> +void object::test<50>() +{ + set_test_name("island in lake"); + + Envelope e(0, 30, 0, 30 ); + Grid ext(e, 10, 10); + + auto g = wkt_reader_.read("MULTIPOLYGON (((5 5, 25 5, 25 25, 5 25, 5 5), (11 11, 11 19, 19 19, 19 11, 11 11)), ((12 12, 14 12, 14 14, 12 14, 12 12)))"); + + auto rci = GridIntersection::getIntersectionFractions(ext, *g); + check_area(*rci, ext, *g); + + auto subd = GridIntersection::subdividePolygon(ext, *g, false); + + check_subdivided_polygon(*g, *subd); +} + +template<> +template<> +void object::test<51>() +{ + set_test_name("subdivide polygon whose edges follow cell boundaries (4)"); + + Envelope e(-180, 180, -90, 90); + Grid ext(e, 0.1, 0.1); + + std::stringstream wkb("0103000000010000004500000039D384ED27C265C03D62F4DC420F3A40E6913F18F8C165C0A702EE79FE143A4020EEEA55E4C165C0BDE2A9471A183A40D6E6FF55C7C165C0AED689CBF11A3A40AF230ED9C0C165C0B491EBA6941B3A40B9E177D3ADC165C087DF4DB7EC1C3A4048A7AE7C96C165C0DF6E490ED81D3A4026E4839E4DC165C0986C3CD8621F3A40F54718062CC165C00000000000203A40B3F0F5B52EC165C00000000000203A40EAEA8EC536C165C00000000000203A40A7936C7539C165C00000000000203A407BF8325104C165C0C0417BF5F1203A40548EC9E2FEC065C00CE544BB0A213A409F3C2CD49AC065C079CBD58F4D223A4037DF88EE59C065C06A17D34CF7223A40FCE07CEA58C065C034A2B437F8223A4000AB23473AC065C096D1C8E715233A40EC18575C1CC065C0C075C58CF0223A402A8E03AF16C065C074EFE192E3223A40CDAB3AAB05C065C0CF4BC5C6BC223A400000000000C065C023BDA8DDAF223A400000000000C065C0FB7953910A233A40CF31207BBDBF65C0395FECBDF8223A401E6CB1DB67BF65C0D40AD3F71A223A40BB99D18F06BF65C0336FD575A8223A40168A743FA7BE65C0ED647094BC223A40ADDEE17668BE65C0A165DD3F16223A40E8137992F4BD65C00000000000203A40D0251C7A0BBE65C00000000000203A4065E3C116BBBD65C062F9F36DC11E3A40CB2F8331A2BD65C0F2E9B12D031E3A40C6A2E9EC64BD65C0C40776FC171C3A4027A25F5B3FBD65C0CF656A12BC193A4035ECF7C43ABD65C04E417E3672193A40F1845E7F12BD65C08E210038F6143A40A585CB2AECBC65C053AEF02E17113A408F52094FE8BC65C0B7ED7BD45F0F3A4030F0DC7B38BD65C0AEB8382A37093A4090A4A487A1BD65C05C1E6B4606053A4043FF04172BBE65C0FE9B1727BE023A409B559FAB2DBE65C0F4893C49BA023A401C261AA4E0BE65C0CEE0EF17B3013A400D37E0F343BF65C00000000000003A40E690D44249BF65C00000000000003A40677DCA3159BF65C00000000000003A4035B6D7825EBF65C00000000000003A405BB395977CBF65C06956B60F79FF3940CAFD0E4581BF65C026C5C72764FF394022A64412BDBF65C0E4D70FB1C1FE3940D0D6C1C1DEBF65C0992842EA76FE39400000000000C065C0FBEAAA402DFE39400000000000C065C0BF61A2410AFE39400000000000C065C0EFACDD76A1FD39400000000000C065C0541B9C887EFD39400952297634C065C0AB96749483FD3940E5F1B4FC40C065C0266DAAEE91FD3940352905DDDEC065C017D68D7747FE39407F8978EBFCC065C097C5C4E6E3FE39406F4BE48233C165C00000000000003A409A7D1EA33CC165C00000000000003A401A14CD0358C165C00000000000003A403925202661C165C00000000000003A4063F2069879C165C0A795422097003A40552FBFD3E4C165C0E46723D74D053A40B8B1D991EAC165C0C7F5EFFACC053A40C26B97361CC265C0F48B12F4170A3A405C1C959B28C265C0329067976F0D3A4039D384ED27C265C03D62F4DC420F3A40"); + auto g = wkb_reader_.readHEX(wkb); + + auto rci = GridIntersection::getIntersectionFractions(ext, *g); + check_area(*rci, ext, *g); + + auto subd = GridIntersection::subdividePolygon(ext, *g, false); + check_subdivided_polygon(*g, *subd); +} + +template<> +template<> +void object::test<52>() +{ + set_test_name("subdivide polygon whose edges follow cell boundaries (5)"); + + Envelope e(-180, 180, -90, 90); + Grid ext(e, 0.05, 0.05); + + std::stringstream wkb("01030000000100000024000000b6b9313d619063c0aeefc34142bc3440306475ab679063c0f12900c633bc3440e695eb6db39063c0b47405db88bb34408f19a88c7f9163c022c7d63384bb3440d591239d019263c06d1b4641f0bc3440e3e13d07169263c07a8ec87729bd3440ec4e779e789263c02ead86c43dbe3440fc51d4997b9263c0897d022846be3440b47405db889263c0265305a392be3440a39410acaa9263c04be658de55bf34401f85eb51b89263c0e1968fa4a4bf3440ef59d768b99263c026732cefaabf34403a1f9e25c89263c00000000000c0344003603c83069363c00000000000c0344044e048a0c19363c00000000000c0344000000000009463c00000000000c0344000000000009463c043723271abc0344000000000009463c0adc090d5adc2344000000000009463c0912a8a5759c334402e1a321e259463c0addba0f65bc3344020b589937b9463c02ff99ffcddc33440ff2268cca49463c04d672783a3c434406bd26d89dc9463c01e1a16a3aec53440151f9f901d9563c0c495b37746c73440761bd47e6b9563c055a52daef1c9344080f44d9a869563c05c74b2d47acb34401213d4f0ad9563c05393e00d69cc3440eeb3ca4ce99563c0527e52edd3cd34406954e0641b9663c0d9ec48f59dcf34402104e44b289663c079060dfd13d0344060f3ead6489663c0ab87f59b29d1344060f3ead6489663c040f1323236e13440172bf28ed68e63c040f1323236e13440172bf28ed68e63c09c4ce66215b8344006ae5057349063c09c4ce66215b83440b6b9313d619063c0aeefc34142bc3440"); + auto g = wkb_reader_.readHEX(wkb); + + auto rci = GridIntersection::getIntersectionFractions(ext, *g); + check_area(*rci, ext, *g); + + auto subd = GridIntersection::subdividePolygon(ext, *g, false); + check_subdivided_polygon(*g, *subd); +} + +template<> +template<> +void object::test<53>() { + set_test_name("includeExterior = true"); + + Envelope e(0, 10, 0, 10); + Grid ext(e, 1, 1); + + auto g = wkt_reader_.read("POLYGON ((9.5 9.5, 11 9.5, 11 11, 9.5 11, 9.5 9.5))"); + + auto subd = GridIntersection::subdividePolygon(ext, *g, true); + + auto inside = wkt_reader_.read("POLYGON ((9.5 9.5, 10 9.5, 10 10, 9.5 10, 9.5 9.5))"); + auto outside = wkt_reader_.read("POLYGON ((9.5 10, 10 10, 10 9.5, 11 9.5, 11 11, 9.5 11, 9.5 10))"); + + ensure_equals(subd->getNumGeometries(), 2u); + + // ensure_equals_geometry would fail here because the generated results have extra nodes at grid cell boundaries + ensure(subd->getGeometryN(0)->equals(inside.get())); + ensure(subd->getGeometryN(1)->equals(outside.get())); +} + } diff --git a/tests/unit/operation/grid/GridTest.cpp b/tests/unit/operation/grid/GridTest.cpp index c8160e3d6..0df5ea8a7 100644 --- a/tests/unit/operation/grid/GridTest.cpp +++ b/tests/unit/operation/grid/GridTest.cpp @@ -283,4 +283,54 @@ void object::test<15>() ensure_equals(g2.getColOffset(g1), 20u); } +template<> +template<> +void object::test<16>() +{ + set_test_name("Crop bounded grid while calculating cell boundaries relative to parent"); + + double res = 1.0/3601; + + Grid grid{global, res, res}; + + const auto env0 = grid.getCellEnvelope(30, 30); + + Envelope cropEnv(global.getMinX() + 17*res + 1e-6, global.getMaxX(), global.getMinY(), global.getMaxY() - 17*res - 1e-6); + + auto cropped = grid.shrinkToFit(cropEnv, false); + + const auto env1 = cropped.getCellEnvelope(13, 13); + + ensure_equals("minX", env0.getMinX(), env1.getMinX()); + ensure_equals("maxX", env0.getMaxX(), env1.getMaxX()); + ensure_equals("minY", env0.getMinY(), env1.getMinY()); + ensure_equals("maxY", env0.getMaxY(), env1.getMaxY()); +} + +template<> +template<> +void object::test<17>() +{ + set_test_name("Crop infinite grid while calculating cell boundaries relative to parent"); + + double res = 1.0/3601; + + Grid grid{global, res, res, global}; + + const auto env0 = grid.getCellEnvelope(30, 30); + + Envelope cropEnv(global.getMinX() + 17*res + 1e-6, global.getMaxX(), global.getMinY(), global.getMaxY() - 17*res - 1e-6); + + auto cropped = grid.shrinkToFit(cropEnv, false); + auto rowOffset = cropped.getRowOffset(grid); + auto colOffset = cropped.getColOffset(grid); + + const auto env1 = cropped.getCellEnvelope(30 - rowOffset, 30 - colOffset); + + ensure_equals("minX", env0.getMinX(), env1.getMinX()); + ensure_equals("maxX", env0.getMaxX(), env1.getMaxX()); + ensure_equals("minY", env0.getMinY(), env1.getMinY()); + ensure_equals("maxY", env0.getMaxY(), env1.getMaxY()); +} + } diff --git a/tests/unit/operation/grid/TraversalAreasTest.cpp b/tests/unit/operation/grid/TraversalAreasTest.cpp index d73cbc3a7..a4664c8e2 100644 --- a/tests/unit/operation/grid/TraversalAreasTest.cpp +++ b/tests/unit/operation/grid/TraversalAreasTest.cpp @@ -2,8 +2,13 @@ #include #include +#include +#include #include +using geos::operation::grid::PerimeterDistance; +using geos::operation::grid::Side; +using geos::operation::grid::Traversal; using geos::operation::grid::TraversalAreas; using geos::geom::CoordinateXY; using geos::geom::Envelope; @@ -11,7 +16,23 @@ using geos::geom::GeometryFactory; namespace tut { struct test_traversalareastest_data : GEOSTestBase { - using TraversalVector = std::vector*>; + using TraversalVector = std::vector; + + Traversal make_traversal(const std::vector& coords, void* parentage=nullptr) + { + Traversal t; + if (!coords.empty()) { + t.enter(coords.front(), Side::NONE, parentage); + } + + for (std::size_t i = 1; i < coords.size() - 1; i++) { + t.add(coords[i]); + } + + t.exit(coords.back(), Side::NONE); + + return t; + } const GeometryFactory& gfact = *GeometryFactory::getDefaultInstance(); }; @@ -30,13 +51,16 @@ void object::test<1>() Envelope b{ 0, 10, 0, 10 }; - std::vector traversal{ { 7, 0 }, { 7, 1 }, { 6, 1 }, { 6, 0 } }; - TraversalVector traversals{ &traversal }; + std::vector coords{ { 7, 0 }, { 7, 1 }, { 6, 1 }, { 6, 0 } }; + Traversal t = make_traversal(coords); + TraversalVector traversals{ &t }; ensure_equals(TraversalAreas::getLeftHandArea(b, traversals), 1); ensure_equals_geometry(TraversalAreas::getLeftHandRings(gfact, b, traversals).get(), "POLYGON ((6 0, 7 0, 7 1, 6 1, 6 0))"); - std::reverse(traversal.begin(), traversal.end()); + std::reverse(coords.begin(), coords.end()); + t = make_traversal(coords); + traversals = { &t }; ensure_equals(TraversalAreas::getLeftHandArea(b, traversals), 99); ensure_equals_geometry(TraversalAreas::getLeftHandRings(gfact, b, traversals).get(), "POLYGON ((0 0, 6 0, 6 1, 7 1, 7 0, 10 0, 10 10, 0 10, 0 0))"); } @@ -49,8 +73,9 @@ void object::test<2>() Envelope b{ 0, 10, 0, 10 }; - std::vector traversal{ { 5, 0 }, { 5, 5 }, { 0, 5 } }; - TraversalVector traversals{ &traversal }; + std::vector coords{ { 5, 0 }, { 5, 5 }, { 0, 5 } }; + Traversal t = make_traversal(coords); + TraversalVector traversals{ &t }; ensure_equals(TraversalAreas::getLeftHandArea(b, traversals), 25); ensure_equals_geometry(TraversalAreas::getLeftHandRings(gfact, b, traversals).get(), "POLYGON ((0 0, 5 0, 5 5, 0 5, 0 0))"); @@ -64,8 +89,9 @@ void object::test<3>() Envelope b{ 0, 10, 0, 10 }; - std::vector traversal{ { 4, 0 }, { 4, 10 } }; - TraversalVector traversals{ &traversal }; + std::vector coords{ { 4, 0 }, { 4, 10 } }; + Traversal t = make_traversal(coords); + TraversalVector traversals{ &t }; ensure_equals(TraversalAreas::getLeftHandArea(b, traversals), 40); ensure_equals_geometry(TraversalAreas::getLeftHandRings(gfact, b, traversals).get(), "POLYGON ((0 0, 4 0, 4 10, 0 10, 0 0))"); @@ -79,9 +105,10 @@ void object::test<4>() Envelope b{ 0, 10, 0, 10 }; - std::vector t1 = { { 2, 10 }, { 2, 0 } }; - std::vector t2 = { { 4, 0 }, { 4, 10 } }; - + std::vector c1 = { { 2, 10 }, { 2, 0 } }; + std::vector c2 = { { 4, 0 }, { 4, 10 } }; + Traversal t1 = make_traversal(c1); + Traversal t2 = make_traversal(c2); TraversalVector traversals{ &t1, &t2 }; ensure_equals(TraversalAreas::getLeftHandArea(b, traversals), 20); @@ -96,21 +123,27 @@ void object::test<5>() Envelope b{ 0, 10, 0, 10 }; - std::vector t1 = { { 2, 0 }, { 2, 2 }, { 0, 2 } }; // 2x2 = 4 - std::vector t2 = { { 3, 10 }, { 3, 0 } }; - std::vector t3 = { { 5, 0 }, { 5, 10 } }; // 2x10 = 20 - std::vector t4 = { { 8, 10 }, { 10, 8 } }; // 2x2/2 = 2 - std::vector t5 = { { 10, 6 }, { 8, 6 }, { 8, 3 }, { 10, 3 } }; // 2x3 = 6 - std::vector t6 = { { 10, 4 }, { 9, 4 }, { 9, 5 }, { 10, 5 } }; // 1x1 = 1 (subtracted) - std::vector t7 = { { 10, 3 }, { 8, 3 }, { 8, 0 } }; // 2x3 = 6 + std::vector c1 = { { 2, 0 }, { 2, 2 }, { 0, 2 } }; // 2x2 = 4 + std::vector c2 = { { 3, 10 }, { 3, 0 } }; + std::vector c3 = { { 5, 0 }, { 5, 10 } }; // 2x10 = 20 + std::vector c4 = { { 8, 10 }, { 10, 8 } }; // 2x2/2 = 2 + std::vector c5 = { { 10, 6 }, { 8, 6 }, { 8, 3 }, { 10, 3 } }; // 2x3 = 6 + std::vector c6 = { { 10, 4 }, { 9, 4 }, { 9, 5 }, { 10, 5 } }; // 1x1 = 1 (subtracted) + std::vector c7 = { { 10, 2 }, { 8, 2 }, { 8, 0 } }; // 2x2 = 4 + + Traversal t1 = make_traversal(c1); + Traversal t2 = make_traversal(c2); + Traversal t3 = make_traversal(c3); + Traversal t4 = make_traversal(c4); + Traversal t5 = make_traversal(c5); + Traversal t6 = make_traversal(c6); + Traversal t7 = make_traversal(c7); TraversalVector traversals{ &t1, &t2, &t3, &t4, &t5, &t6, &t7 }; - ensure_equals(TraversalAreas::getLeftHandArea(b, traversals), 4 + 20 + 2 + 6 - 1 + 6); - ensure_equals_geometry(TraversalAreas::getLeftHandRings(gfact, b, traversals).get(), "MULTIPOLYGON (((2 0, 2 2, 0 2, 0 0, 2 0)), " - "((3 10, 3 0, 5 0, 5 10, 3 10)), " - "((8 10, 10 8, 10 10, 8 10)), " - "((10 6, 8 6, 8 3, 10 3, 10 3, 8 3, 8 0, 10 0, 10 4, 9 4, 9 5, 10 5, 10 6)))"); + ensure_equals(TraversalAreas::getLeftHandArea(b, traversals), 4 + 20 + 2 + 6 - 1 + 4); + + ensure_equals_geometry(TraversalAreas::getLeftHandRings(gfact, b, traversals).get(), "MULTIPOLYGON (((2 0, 2 2, 0 2, 0 0, 2 0)), ((3 10, 3 0, 5 0, 5 10, 3 10)), ((8 10, 10 8, 10 10, 8 10)), ((10 6, 8 6, 8 3, 10 3, 10 4, 9 4, 9 5, 10 5, 10 6)), ((10 2, 8 2, 8 0, 10 0, 10 2)))"); } template<> @@ -135,7 +168,8 @@ void object::test<7>() Envelope b{ 0, 10, 0, 10 }; - std::vector t1{ { 4, 0 }, { 4, 0 } }; + std::vector c1{ { 4, 0 }, { 4, 0 } }; + Traversal t1 = make_traversal(c1); TraversalVector traversals{ &t1 }; ensure_THROW(TraversalAreas::getLeftHandArea(b, traversals), std::exception); @@ -150,7 +184,8 @@ void object::test<8>() Envelope b{ 0, 10, 0, 10 }; - std::vector t1 = { { 1, 1 }, { 2, 1 }, { 2, 2 }, { 1, 2 }, { 1, 1 } }; + std::vector c1 = { { 1, 1 }, { 2, 1 }, { 2, 2 }, { 1, 2 }, { 1, 1 } }; + Traversal t1 = make_traversal(c1); TraversalVector traversals{ &t1 }; ensure_equals(TraversalAreas::getLeftHandArea(b, traversals), 1); @@ -161,11 +196,12 @@ template<> template<> void object::test<9>() { - set_test_name("Closed ring ccw overlapping edge"); + set_test_name("Closed ring ccw touching edge"); Envelope b{ 0, 10, 00, 10 }; - std::vector t1 = { { 1, 0 }, { 2, 1 }, { 1, 1 }, { 1, 0 } }; + std::vector c1 = { { 1, 0 }, { 2, 1 }, { 1, 1 }, { 1, 0 } }; + Traversal t1 = make_traversal(c1); TraversalVector traversals{ &t1 }; ensure_equals(TraversalAreas::getLeftHandArea(b, traversals), 0.5); @@ -180,10 +216,12 @@ void object::test<10>() Envelope b{ 0, 10, 0, 10 }; - std::vector t1 = { { 1, 1 }, { 1, 2 }, { 2, 2 }, { 2, 1 }, { 1, 1 } }; + std::vector c1 = { { 1, 1 }, { 1, 2 }, { 2, 2 }, { 2, 1 }, { 1, 1 } }; + Traversal t1 = make_traversal(c1); TraversalVector traversals{ &t1 }; ensure_equals(TraversalAreas::getLeftHandArea(b, traversals), 99); + ensure_equals_geometry(TraversalAreas::getLeftHandRings(gfact, b, traversals).get(), "POLYGON ((0 0, 10 0, 10 10, 0 10, 0 0), (1 1, 1 2, 2 2, 2 1, 1 1))"); } @@ -195,12 +233,14 @@ void object::test<11>() Envelope b{ 0, 10, 0, 10 }; - std::vector t1 = { { 1, 1 }, { 1, 2 }, { 2, 2 }, { 2, 1 }, { 1, 1 } }; - std::vector t2 = { { 10, 5 }, { 10, 5 } }; + std::vector c1 = { { 1, 1 }, { 1, 2 }, { 2, 2 }, { 2, 1 }, { 1, 1 } }; + std::vector c2 = { { 10, 5 }, { 10, 5 } }; + Traversal t1 = make_traversal(c1); + Traversal t2 = make_traversal(c2); TraversalVector traversals{ &t1, &t2 }; ensure_equals(TraversalAreas::getLeftHandArea(b, traversals), 99); - ensure_equals_geometry(TraversalAreas::getLeftHandRings(gfact, b, traversals).get(), "POLYGON ((0 0, 10 0, 10 10, 0 10, 0 0), (1 1, 1 2, 2 2, 2 1, 1 1))"); + ensure_equals_geometry(TraversalAreas::getLeftHandRings(gfact, b, traversals).get(), "POLYGON ((0 0, 10 0, 10 5, 10 10, 0 10, 0 0), (1 1, 1 2, 2 2, 2 1, 1 1))"); } template<> @@ -211,7 +251,8 @@ void object::test<12>() Envelope b{ 0, 10, 0, 10 }; - std::vector t1 = { { 0, 0 }, { 2, 2 }, { 3, 2 }, { 0, 0 } }; + std::vector c1 = { { 0, 0 }, { 2, 2 }, { 3, 2 }, { 0, 0 } }; + Traversal t1 = make_traversal(c1); TraversalVector traversals{ &t1 }; ensure_equals(TraversalAreas::getLeftHandArea(b, traversals), 99); @@ -224,13 +265,14 @@ void object::test<13>() { set_test_name("Closed ring cw touching edge interior"); - Envelope b{ 0, 10, 0, 10 }; - - std::vector t1 = { { 1, 0 }, { 2, 2 }, { 3, 2 }, { 1, 0 } }; + Envelope b(0, 10, 0, 10); + std::vector c1{{8, 8}, {5, 0}, {5, 8}, {8, 8}}; + Traversal t1 = make_traversal(c1); TraversalVector traversals{ &t1 }; - ensure_equals(TraversalAreas::getLeftHandArea(b, traversals), 99); - ensure_equals_geometry(TraversalAreas::getLeftHandRings(gfact, b, traversals).get(), "POLYGON ((0 0, 10 0, 10 10, 0 10, 0 0), (1 0, 2 2, 3 2, 1 0))"); + ensure_equals(TraversalAreas::getLeftHandArea(b, traversals), 88); + ensure_equals_geometry(TraversalAreas::getLeftHandRings(gfact, b, traversals).get(), + "POLYGON ((0 10, 10 10, 10 0, 5 0, 0 0, 0 10), (8 8, 5 8, 5 0, 8 8))"); } template<> @@ -239,13 +281,14 @@ void object::test<14>() { set_test_name("Closed ring cw overlapping edge"); - Envelope b{ 0, 10, 0, 10 }; - - std::vector t1 = { { 1, 0 }, { 1, 1 }, { 2, 1 }, { 1, 0 } }; + Envelope b(0, 10, 0, 10); + std::vector c1{{8, 8}, {8, 0}, {5, 0}, {5, 8}, {8, 8}}; + Traversal t1 = make_traversal(c1); TraversalVector traversals{ &t1 }; - ensure_equals(TraversalAreas::getLeftHandArea(b, traversals), 99.5); - ensure_equals_geometry(TraversalAreas::getLeftHandRings(gfact, b, traversals).get(), "POLYGON ((0 0, 10 0, 10 10, 0 10, 0 0), (1 0, 1 1, 2 1, 1 0))"); + ensure_equals(TraversalAreas::getLeftHandArea(b, traversals), 76); + ensure_equals_geometry(TraversalAreas::getLeftHandRings(gfact, b, traversals).get(), + "POLYGON ((0 10, 10 10, 10 0, 8 0, 8 8, 5 8, 5 0, 0 0, 0 10))"); } template<> @@ -256,12 +299,12 @@ void object::test<15>() Envelope b{ 0, 10, 0, 10 }; - std::vector t1{ { 4, 0 }, { 10, 0 } }; + std::vector c1{ { 4, 0 }, { 10, 0 } }; + Traversal t1 = make_traversal(c1); TraversalVector traversals{ &t1 }; ensure_equals(TraversalAreas::getLeftHandArea(b, traversals), 100); - // FIXME remove repeated point from result - ensure_equals_geometry(TraversalAreas::getLeftHandRings(gfact, b, traversals).get(), "POLYGON ((0 0, 4 0, 10 0, 10 0, 10 10, 0 10, 0 0))"); + ensure_equals_geometry(TraversalAreas::getLeftHandRings(gfact, b, traversals).get(), "POLYGON ((0 0, 4 0, 10 0, 10 10, 0 10, 0 0))"); } template<> @@ -272,7 +315,8 @@ void object::test<16>() Envelope b{ 2, 3, 2, 3 }; - std::vector t1{ { 2, 2 }, { 2, 2.5 }, { 2, 2.5 } }; + std::vector c1{ { 2, 2 }, { 2, 2.5 }, { 2, 2.5 } }; + Traversal t1 = make_traversal(c1); TraversalVector traversals{ &t1 }; ensure_equals(TraversalAreas::getLeftHandArea(b, traversals), 0); @@ -281,18 +325,281 @@ void object::test<16>() template<> template<> -void object::test<17>() { +void object::test<17>() +{ + set_test_name("interior and edge traversal"); + + Envelope b(6, 7, 3, 4); + std::vector c1{{7, 3}, {6, 4}, {7, 4}}; + Traversal t1 = make_traversal(c1); + TraversalVector traversals{ &t1 }; + + ensure_equals(TraversalAreas::getLeftHandArea(b, traversals), 0.5); + ensure_equals_geometry(TraversalAreas::getLeftHandRings(gfact, b, traversals).get(), "POLYGON ((6 3, 7 3, 6 4, 6 3))"); +} + +template<> +template<> +void object::test<18>() +{ + set_test_name("interior-edge segment-interior traversal"); + + Envelope b(0, 10, 0, 10); + std::vector c1{{10, 5}, {8, 0}, {4, 0}, {0, 3}}; + Traversal t1 = make_traversal(c1); + TraversalVector traversals{ &t1 }; + + ensure_equals(TraversalAreas::getLeftHandArea(b, traversals), 11); + ensure_equals_geometry(TraversalAreas::getLeftHandRings(gfact, b, traversals).get(), "MULTIPOLYGON (((10 5, 8 0, 10 0, 10 5)), ((4 0, 0 3, 0 0, 4 0)))"); +} + +template<> +template<> +void object::test<19>() +{ + set_test_name("interior-edge point-interior traversal"); + + Envelope b(0, 10, 0, 10); + std::vector c1{{10, 5}, {8, 0}, {0, 3}}; + Traversal t1 = make_traversal(c1); + TraversalVector traversals{ &t1 }; + + ensure_equals(TraversalAreas::getLeftHandArea(b, traversals), 17); + ensure_equals_geometry(TraversalAreas::getLeftHandRings(gfact, b, traversals).get(), "MULTIPOLYGON (((8 0, 10 5, 10 0, 8 0)), ((8 0, 0 0, 0 3, 8 0)))"); +} + +template<> +template<> +void object::test<20>() +{ + set_test_name("interior-edge point-interior traversal, with repeated points"); + + Envelope b(0, 10, 0, 10); + std::vector c1{{10, 5}, {8, 0}, {8, 0}, {0, 3}}; + Traversal t1 = make_traversal(c1); + TraversalVector traversals{ &t1 }; + + ensure_equals(TraversalAreas::getLeftHandArea(b, traversals), 17); + ensure_equals_geometry(TraversalAreas::getLeftHandRings(gfact, b, traversals).get(), "MULTIPOLYGON (((8 0, 10 5, 10 0, 8 0)), ((8 0, 0 0, 0 3, 8 0)))"); +} + +template<> +template<> +void object::test<21>() +{ + set_test_name("bouncing off multiple edges"); + + Envelope b(0, 10, 0, 10); + + std::vector c1{{10, 5}, {5, 10}, {2, 0}, {0, 5}}; + Traversal t1 = make_traversal(c1); + TraversalVector traversals{ &t1 }; + + ensure_equals(TraversalAreas::getLeftHandArea(b, traversals), 57.5); + ensure_equals_geometry(TraversalAreas::getLeftHandRings(gfact, b, traversals).get(), "MULTIPOLYGON (((10 5, 5 10, 2 0, 10 0, 10 5)), ((2 0, 0 5, 0 0, 2 0)))"); +} + +template<> +template<> +void object::test<22>() +{ + set_test_name("enter top, bounce bottom, exit right"); + + Envelope b(0, 10, 0, 10); + + std::vector c1{{5, 10}, {5, 0}, {10, 2}}; + Traversal t1 = make_traversal(c1); + TraversalVector traversals{ &t1 }; + + ensure_equals(TraversalAreas::getLeftHandArea(b, traversals), 45); + ensure_equals_geometry(TraversalAreas::getLeftHandRings(gfact, b, traversals).get(), + "POLYGON ((5 10, 5 0, 10 2, 10 10, 5 10))" + ); +} + +template<> +template<> +void object::test<23>() +{ + set_test_name("along top, then bottom to top"); + + Envelope b(0, 10, 0, 10); + std::vector c1{{5, 10}, {0, 10}}; + std::vector c2{{5, 0}, {5, 10}}; + Traversal t1 = make_traversal(c1); + Traversal t2 = make_traversal(c2); + TraversalVector traversals{ &t1, &t2 }; + + ensure_equals(TraversalAreas::getLeftHandArea(b, traversals), 50); + ensure_equals_geometry(TraversalAreas::getLeftHandRings(gfact, b, traversals).get(), + "POLYGON ((5 10, 0 10, 0 0, 5 0, 5 10))"); +} + +template<> +template<> +void object::test<24>() +{ + set_test_name("two traversals, touching in interior"); + // This would occur when two touching holes are present in the cell, and the shell is not present. + // Or it could occur when a shell and hole touch + + Envelope b(0, 10, 0, 10); + std::vector c1{{10, 0}, {5, 5}, {10, 10}}; + std::vector c2{{0, 10}, {5, 5}, {0, 0}}; + Traversal t1 = make_traversal(c1, &c1); + Traversal t2 = make_traversal(c2, &c2); + TraversalVector traversals{ &t1, &t2 }; + + ensure_equals(TraversalAreas::getLeftHandArea(b, traversals), 50); + ensure_equals_geometry(TraversalAreas::getLeftHandRings(gfact, b, traversals).get(), + "MULTIPOLYGON (((5 5, 0 10, 10 10, 5 5)), ((5 5, 10 0, 0 0, 5 5)))"); +} + +template<> +template<> +void object::test<25>() +{ + set_test_name("two traversals, touching in interior, plus complete hole"); + // This is the same as #24 but the complete hole forces use of the polygonizer + + Envelope b(0, 10, 0, 10); + std::vector c1{{10, 0}, {5, 5}, {10, 10}}; + std::vector c2{{0, 10}, {5, 5}, {0, 0}}; + std::vector c3{{5, 1}, {5, 2}, {6, 2}, {6, 1}, {5, 1}}; + Traversal t1 = make_traversal(c1, &c1); + Traversal t2 = make_traversal(c2, &c2); + Traversal t3 = make_traversal(c3, &c3); + TraversalVector traversals{ &t1, &t2, &t3 }; + + ensure_equals(TraversalAreas::getLeftHandArea(b, traversals), 49); + ensure_equals_geometry(TraversalAreas::getLeftHandRings(gfact, b, traversals).get(), + "MULTIPOLYGON (((0 10, 10 10, 5 5, 0 10)), ((0 0, 5 5, 10 0, 0 0), (5 1, 6 1, 6 2, 5 2, 5 1)))"); +} + +template<> +template<> +void object::test<26>() +{ + set_test_name("multiple holes touching at endpoints"); + + Envelope b(0, 10, 0, 10); + std::vector c1{{0, 10}, {3, 1}, {0, 0}}; + std::vector c2{{0, 10}, {10, 9}}; + std::vector c3{{10, 7}, {0, 10}}; + Traversal t1 = make_traversal(c1, &c1); + Traversal t2 = make_traversal(c2, &c2); + Traversal t3 = make_traversal(c3, &c2); + TraversalVector traversals{ &t1, &t2, &t3 }; + + ensure_equals(TraversalAreas::getLeftHandArea(b, traversals), 75); + ensure_equals_geometry(TraversalAreas::getLeftHandRings(gfact, b, traversals).get(), + "MULTIPOLYGON (((0 10, 3 1, 0 0, 10 0, 10 7, 0 10)), ((0 10, 10 9, 10 10, 0 10)))"); +} + +template<> +template<> +void object::test<27>() +{ + set_test_name("island with lake"); + + Envelope b(0, 10, 0, 10); + std::vector c1{{1, 1}, {9, 1}, {9, 9}, {1, 9}, {1, 1}}; + std::vector c2{{2, 2}, {2, 4}, {4, 4}, {4, 2}, {2, 2}}; + Traversal island = make_traversal(c1, &c1); + Traversal lake = make_traversal(c2, &c2); + TraversalVector traversals{ &island, &lake }; + + ensure_equals(TraversalAreas::getLeftHandArea(b, traversals), 60); + ensure_equals_geometry(TraversalAreas::getLeftHandRings(gfact, b, traversals).get(), + "POLYGON ((1 1, 1 9, 9 9, 9 1, 1 1), (2 2, 4 2, 4 4, 2 4, 2 2))"); +} + +template<> +template<> +void object::test<28>() +{ set_test_name("lake with island"); Envelope b(0, 10, 0, 10); - std::vector lake{{1, 1}, {1, 9}, {9, 9}, {9, 1}, {1, 1}}; - std::vector island{{2, 2}, {4, 2}, {4, 4}, {2, 4}, {2, 2}}; + std::vector c1{{1, 1}, {1, 9}, {9, 9}, {9, 1}, {1, 1}}; + std::vector c2{{2, 2}, {4, 2}, {4, 4}, {2, 4}, {2, 2}}; + Traversal lake = make_traversal(c1, &c1); + Traversal island = make_traversal(c2, &c2); TraversalVector traversals{ &island, &lake }; ensure_equals(TraversalAreas::getLeftHandArea(b, traversals), 40); - //ensure_equals_geometry(TraversalAreas::getLeftHandRings(gfact, b, traversals).get(), - //"MULTIPOLYGON (((0 0, 10 0, 10 10, 0 10, 0 0), (1 1, 9 1, 9 9, 1 9, 1 1)), ((2 2, 4 2, 4 4, 2 4, 2 2)))"); + ensure_equals_geometry(TraversalAreas::getLeftHandRings(gfact, b, traversals).get(), + "MULTIPOLYGON (((0 0, 10 0, 10 10, 0 10, 0 0), (1 1, 9 1, 9 9, 1 9, 1 1)), ((2 2, 4 2, 4 4, 2 4, 2 2)))"); } +template<> +template<> +void object::test<29>() +{ + set_test_name("traversals define region with no area"); + // this can occur only with invalid polygons but we have historically handled these correctly + + Envelope b(0, 10, 0, 10); + std::vector c1{{0, 5}, {10, 5}}; + std::vector c2{{10, 5}, {0, 5}}; + Traversal t1 = make_traversal(c1, &c1); + Traversal t2 = make_traversal(c2, &c2); + TraversalVector traversals{ &t1, &t2 }; + + ensure_equals(TraversalAreas::getLeftHandArea(b, traversals), 0); +} + +template<> +template<> +void object::test<30>() +{ + set_test_name("PerimeterDistance::isLessThan robustness test"); + + Envelope env{ 146.52500000000001, 146.57500000000002, -33.625, -33.575000000000003}; + + { + // ----p1--p2--+ + // | + // | + CoordinateXY p1{std::nextafter(env.getMaxX(), geos::DoubleNegInfinity), env.getMaxY()}; + CoordinateXY p2{env.getMaxX(), env.getMaxY()}; + + ensure("points along ymax", PerimeterDistance::isLessThan(env, p1, p2)); + } + + { + // | + // | + // ----p2--p1--+ + CoordinateXY p1{env.getMaxX(), env.getMinY()}; + CoordinateXY p2{std::nextafter(env.getMaxX(), geos::DoubleNegInfinity), env.getMinY()}; + + ensure("points along ymin", PerimeterDistance::isLessThan(env, p1, p2)); + } + + { + // +------------ + // p2 + // | + // p1 + // | + CoordinateXY p1{env.getMinX(), std::nextafter(env.getMaxY(), geos::DoubleNegInfinity)}; + CoordinateXY p2{env.getMinX(), env.getMaxY()}; + + ensure("points along xmin", PerimeterDistance::isLessThan(env, p1, p2)); + } + + { + // p1 + // | + // p2 + // | + // -----------+ + CoordinateXY p1{env.getMaxX(), std::nextafter(env.getMinY(), geos::DoubleInfinity)}; + CoordinateXY p2{env.getMaxX(), env.getMinY()}; + + ensure("points along xmax", PerimeterDistance::isLessThan(env, p1, p2)); + } +} } diff --git a/tests/unit/utility.h b/tests/unit/utility.h index 5d3993d10..562a95dc8 100644 --- a/tests/unit/utility.h +++ b/tests/unit/utility.h @@ -581,6 +581,14 @@ struct wkb_hex_decoder { // Base class with helpers struct GEOSTestBase { + static void + skip_on_ci() + { + if (std::getenv("CI")) { + skip("Skip slow test on CI"); + } + } + static std::string load_resource(const std::string& fname) { #if HAVE_STD_FILESYSTEM if (RESOURCE_DIR == "") {