diff --git a/src/spatial/modules/geos/geos_geometry.hpp b/src/spatial/modules/geos/geos_geometry.hpp index a275e915..36d66961 100644 --- a/src/spatial/modules/geos/geos_geometry.hpp +++ b/src/spatial/modules/geos/geos_geometry.hpp @@ -15,6 +15,7 @@ class GeosGeometry { public: // constructor GeosGeometry(GEOSContextHandle_t handle_p, GEOSGeometry *geom_p); + GeosGeometry(GEOSContextHandle_t handle_p, double xmin, double ymin, double xmax, double ymax); // disable copy GeosGeometry(const GeosGeometry &) = delete; @@ -40,6 +41,7 @@ class GeosGeometry { bool is_valid() const; bool is_empty() const; + GeosGeometry get_clone() const; GeosGeometry get_boundary() const; GeosGeometry get_centroid() const; GeosGeometry get_convex_hull() const; @@ -89,6 +91,7 @@ class GeosGeometry { // default tolerance is max(height/width) / 1000 GeosGeometry get_maximum_inscribed_circle() const; GeosGeometry get_point_n(int n) const; + GeosGeometry get_geometry_n(int n) const; GeosGeometry get_linemerged(bool directed) const; GeosGeometry get_concave_hull(const double ratio, const bool allowHoles) const; @@ -105,6 +108,13 @@ class GeosGeometry { void get_extent(double &xmin, double &ymin, double &xmax, double &ymax) const { GEOSGeom_getExtent_r(handle, geom, &xmin, &ymin, &xmax, &ymax); } + size_t get_dimension() const { + return GEOSGeom_getDimensions_r(handle, geom); + } + size_t get_num_geometries() const { + return GEOSGetNumGeometries_r(handle, geom); + } + size_t get_num_vertices() const; private: GEOSContextHandle_t handle; @@ -208,6 +218,9 @@ class PreparedGeosGeometry { //-- GeosGeometry --// inline GeosGeometry::GeosGeometry(GEOSContextHandle_t handle_p, GEOSGeometry *geom_p) : handle(handle_p), geom(geom_p) { } +inline GeosGeometry::GeosGeometry(GEOSContextHandle_t handle_p, double xmin, double ymin, double xmax, double ymax) + : handle(handle_p), geom(GEOSGeom_createRectangle_r(handle_p, xmin, ymin, xmax, ymax)) { +} inline GeosGeometry::GeosGeometry(GeosGeometry &&other) noexcept : handle(other.handle), geom(other.geom) { other.geom = nullptr; other.handle = nullptr; @@ -296,6 +309,10 @@ inline bool GeosGeometry::is_empty() const { return GEOSisEmpty_r(handle, geom); } +inline GeosGeometry GeosGeometry::get_clone() const { + return GeosGeometry(handle, GEOSGeom_clone_r(handle, geom)); +} + inline GeosGeometry GeosGeometry::get_boundary() const { return GeosGeometry(handle, GEOSBoundary_r(handle, geom)); } @@ -397,6 +414,13 @@ inline GeosGeometry GeosGeometry::get_point_n(int n) const { return GeosGeometry(handle, point); } +inline GeosGeometry GeosGeometry::get_geometry_n(int n) const { + // TODO: GEOSGeometryN_r returns a pointer into the internal storage + // of geom, so we need to return a copy if we want to keep the interface of GeosGeometry the same + auto subgeom = GEOSGetGeometryN_r(handle, geom, n); + return GeosGeometry(handle, GEOSGeom_clone_r(handle, subgeom)); +} + inline bool GeosGeometry::contains(const GeosGeometry &other) const { return GEOSContains_r(handle, geom, other.geom); } @@ -524,6 +548,69 @@ inline PreparedGeosGeometry GeosGeometry::get_prepared() const { return PreparedGeosGeometry(handle, *this); } +inline size_t GeosGeometry::get_num_vertices() const { + size_t num_vertices = 0; + + if (GEOSisEmpty_r(handle, geom)) { + return num_vertices; + } + + switch (GEOSGeomTypeId_r(handle, geom)) { + case GEOS_POINT: { + num_vertices = 1; + break; + }; + case GEOS_LINESTRING: + case GEOS_LINEARRING: { + const int line_vertecies = GEOSGeomGetNumPoints_r(handle, geom); + if (line_vertecies == -1) { + throw InvalidInputException("Could not get number of points for LineString or LinearRing input"); + } + + num_vertices = static_cast(line_vertecies); + break; + } + case GEOS_POLYGON: { + const int exterior_ring_vertecies = GEOSGeomGetNumPoints_r(handle, GEOSGetExteriorRing_r(handle, geom)); + if (exterior_ring_vertecies == -1) { + throw InvalidInputException("Could not get number of points for polygon exterior ring"); + } + num_vertices += static_cast(exterior_ring_vertecies); + + for (size_t i = 0; i < GEOSGetNumInteriorRings_r(handle, geom); i++) { + const int interior_ring_vertecies = GEOSGeomGetNumPoints_r(handle, GEOSGetInteriorRingN_r(handle, geom, i)); + if (interior_ring_vertecies == -1) { + throw InvalidInputException("Could not get number of points for polygon interior ring %d", i); + } + num_vertices += static_cast(interior_ring_vertecies); + } + break; + } + case GEOS_MULTIPOINT: { + // For a MultiPoint, no need to check nested geometries as there are none.. + num_vertices = GEOSGetNumGeometries_r(handle, geom); + break; + } + case GEOS_MULTILINESTRING: + case GEOS_MULTIPOLYGON: + case GEOS_GEOMETRYCOLLECTION: { + for (size_t i = 0; i < GEOSGetNumGeometries_r(handle, geom); i++) { + // TODO: Ideally, there should be a read-only reference into the geometry that we receive + // instead of cloning it.. + GeosGeometry subgeom {handle, GEOSGeom_clone_r(handle, GEOSGetGeometryN_r(handle, geom, i))}; + num_vertices += subgeom.get_num_vertices(); + } + break; + } + default: { + throw InvalidInputException("Geometry type not implemented for get_num_vertices: %d", + GEOSGeomTypeId_r(handle, geom)); + } + } + + return num_vertices; +} + //-- PreparedGeosGeometry --// inline bool PreparedGeosGeometry::contains(const GeosGeometry &other) const { diff --git a/src/spatial/modules/geos/geos_module.cpp b/src/spatial/modules/geos/geos_module.cpp index ee75015a..5350f5ef 100644 --- a/src/spatial/modules/geos/geos_module.cpp +++ b/src/spatial/modules/geos/geos_module.cpp @@ -2096,6 +2096,153 @@ struct ST_SimplifyPreserveTopology { } }; +struct ST_Subdivide { + static void SubdivideRecursive(GEOSContextHandle_t handle_p, GeosCollection &collection, const GeosGeometry &geom, + const size_t max_vertices, const size_t dimension, const size_t depth) { + const size_t max_depth = 50; + + // If we get a lower-dimensional object from an intersection, we abort + if (geom.get_dimension() < dimension) { + return; + } + + // A MultiPoint is ignored here on purpose as MultiPoints get treated as one + // object compared to multiple distinct objects + if (geom.type() >= GEOS_MULTILINESTRING && geom.type() <= GEOS_GEOMETRYCOLLECTION) { + for (size_t i = 0; i < geom.get_num_geometries(); i++) { + const auto subgeom = geom.get_geometry_n(i); + + // Do not increment depth as we are still on the same level, just processing individual + // parts of the geometry + SubdivideRecursive(handle_p, collection, subgeom, max_vertices, dimension, depth); + } + return; + } + + if (geom.get_num_vertices() <= max_vertices) { + collection.add(std::move(geom.get_clone())); + return; + } + + // Went so far that we will just add the rest all at once. + if (depth > max_depth) { + collection.add(std::move(geom.get_clone())); + return; + } + + double xmin, ymin, xmax, ymax; + geom.get_extent(xmin, ymin, xmax, ymax); + + const double width = xmax - xmin; + const double height = ymax - ymin; + + if (width == 0.0 && height == 0.0) { + if (geom.type() == GEOS_POINT) { + collection.add(std::move(geom.get_clone())); + } + + return; + } + + if (width == 0.0) { + xmin -= std::numeric_limits::max(); + xmax += std::numeric_limits::max(); + } + + if (height == 0.0) { + ymin -= std::numeric_limits::max(); + ymax += std::numeric_limits::max(); + } + + double xmin_a, xmax_a, ymin_a, ymax_a; + double xmin_b, xmax_b, ymin_b, ymax_b; + if (width > height) { + xmin_a = xmin; + xmax_a = (xmax + xmin) / 2.0; + ymin_a = ymin; + ymax_a = ymax; + + xmin_b = (xmax + xmin) / 2.0; + xmax_b = xmax; + ymin_b = ymin; + ymax_b = ymax; + } else { + xmin_a = xmin; + xmax_a = xmax; + ymin_a = ymin; + ymax_a = (ymax + ymin) / 2.0; + + xmin_b = xmin; + xmax_b = xmax; + ymin_b = (ymax + ymin) / 2.0; + ymax_b = ymax; + } + + { + const GeosGeometry clipping_rect_a = GeosGeometry(handle_p, xmin_a, ymin_a, xmax_a, ymax_a); + const GeosGeometry clipped_a = geom.get_intersection(clipping_rect_a); + if (!clipped_a.is_empty()) { + SubdivideRecursive(handle_p, collection, clipped_a, max_vertices, dimension, depth + 1); + } + } + { + const GeosGeometry clipping_rect_b = GeosGeometry(handle_p, xmin_b, ymin_b, xmax_b, ymax_b); + const GeosGeometry clipped_b = geom.get_intersection(clipping_rect_b); + if (!clipped_b.is_empty()) { + SubdivideRecursive(handle_p, collection, clipped_b, max_vertices, dimension, depth + 1); + } + } + + return; + } + + static void Execute(DataChunk &args, ExpressionState &state, Vector &result) { + auto &lstate = LocalState::ResetAndGet(state); + + BinaryExecutor::Execute( + args.data[0], args.data[1], result, args.size(), + [&](const string_t &geom_blob, const uint32_t max_vertices) { + if (max_vertices < 5) { + throw InvalidInputException("max_vertices needs to be larger or equal to 5"); + } + + const auto geom = lstate.Deserialize(geom_blob); + + if (geom.type() == GEOS_GEOMETRYCOLLECTION) { + throw InvalidInputException("Cannot subdivide GeometryCollection"); + } + + if (geom.is_empty()) { + return lstate.Serialize(result, geom); + } + + GeosCollection collection {lstate.GetContext()}; + SubdivideRecursive(lstate.GetContext(), collection, geom, max_vertices, geom.get_dimension(), 0); + + return lstate.Serialize(result, collection.get_collection()); + }); + } + + static void Register(ExtensionLoader &loader) { + FunctionBuilder::RegisterScalar(loader, "ST_Subdivide", [](ScalarFunctionBuilder &func) { + func.AddVariant([](ScalarFunctionVariantBuilder &variant) { + variant.AddParameter("geom", LogicalType::GEOMETRY()); + variant.AddParameter("max_vertices", LogicalType::UINTEGER); + variant.SetReturnType(LogicalType::GEOMETRY()); + + variant.SetInit(LocalState::Init); + variant.SetFunction(Execute); + }); + + func.SetDescription( + "Recursively splits a geometry into sub-geometries until the number of vertices of each are below the " + "threshold given by max_vertices. Accepts any type of input except for a GeometryCollection."); + func.SetTag("ext", "spatial"); + func.SetTag("category", "relation"); + }); + } +}; + struct ST_Touches : SymmetricPreparedBinaryFunction { static bool ExecutePredicateNormal(const GeosGeometry &lhs, const GeosGeometry &rhs) { return lhs.touches(rhs); @@ -3029,6 +3176,7 @@ void RegisterGEOSModule(ExtensionLoader &loader) { ST_ShortestLine::Register(loader); ST_Simplify::Register(loader); ST_SimplifyPreserveTopology::Register(loader); + ST_Subdivide::Register(loader); ST_Touches::Register(loader); ST_Union::Register(loader); ST_VoronoiDiagram::Register(loader); diff --git a/test/sql/geos/st_subdivide.test b/test/sql/geos/st_subdivide.test new file mode 100644 index 00000000..eb2407a7 --- /dev/null +++ b/test/sql/geos/st_subdivide.test @@ -0,0 +1,41 @@ +require spatial + +query I +SELECT st_astext(st_subdivide(st_geomfromtext('POINT EMPTY'), 5)); +---- +POINT EMPTY + +query I +SELECT st_astext(st_subdivide(st_geomfromtext('POINT (1 1)'), 5)); +---- +GEOMETRYCOLLECTION (POINT (1 1)) + +query I +select st_subdivide(st_geomfromtext('MULTIPOINT ((1 1), (2 2), (3 3), (4 4), (5 5), (6 6))'), 5); +---- +GEOMETRYCOLLECTION (MULTIPOINT (1 1, 2 2, 3 3), MULTIPOINT (4 4, 5 5, 6 6)) + +query I +SELECT st_astext(st_subdivide(st_geomfromtext('LINESTRING (1 1, 2 2, 3 3, 4 4, 5 5, 6 6, 7 7, 8 8)'), 5)); +---- +GEOMETRYCOLLECTION (LINESTRING (1 1, 2 2, 3 3, 4 4, 4.5 4.5), LINESTRING (4.5 4.5, 5 5, 6 6, 7 7, 8 8)) + +query I +select st_astext(st_subdivide(st_geomfromtext('MULTILINESTRING ((1 1, 2 2, 3 3, 4 4, 5 5, 6 6), (7 7, 8 8, 9 9, 10 10, 11 11, 12 12, 13 13))'), 5)); +---- +GEOMETRYCOLLECTION (LINESTRING (1 1, 2 2, 3 3, 3.5 3.5), LINESTRING (3.5 3.5, 4 4, 5 5, 6 6), LINESTRING (7 7, 8 8, 9 9, 10 10), LINESTRING (10 10, 11 11, 12 12, 13 13)) + +query I +SELECT st_astext(st_subdivide(st_geomfromtext('POLYGON((0 0, 0 10, 5 8, 10 10, 10 0, 5 2, 0 0))'), 7)); +---- +GEOMETRYCOLLECTION (POLYGON ((0 0, 0 10, 5 8, 10 10, 10 0, 5 2, 0 0))) + +query I +SELECT st_astext(st_subdivide(st_geomfromtext('POLYGON((0 0, 0 10, 5 8, 10 10, 10 0, 5 2, 0 0))'), 5)); +---- +GEOMETRYCOLLECTION (POLYGON ((5 5, 5 2, 0 0, 0 5, 5 5)), POLYGON ((10 5, 10 0, 5 2, 5 5, 10 5)), POLYGON ((5 8, 5 5, 0 5, 0 10, 5 8)), POLYGON ((10 10, 10 5, 5 5, 5 8, 10 10))) + +query I +select st_astext(st_subdivide(st_geomfromtext('MULTIPOLYGON(((0 0, 0 10, 5 8, 10 10, 10 0, 5 2, 0 0)), ((20 20, 20 30, 25 28, 30 30, 30 20, 25 22, 20 20)))'), 5)); +---- +GEOMETRYCOLLECTION (POLYGON ((5 5, 5 2, 0 0, 0 5, 5 5)), POLYGON ((10 5, 10 0, 5 2, 5 5, 10 5)), POLYGON ((5 8, 5 5, 0 5, 0 10, 5 8)), POLYGON ((10 10, 10 5, 5 5, 5 8, 10 10)), POLYGON ((25 25, 25 22, 20 20, 20 25, 25 25)), POLYGON ((30 25, 30 20, 25 22, 25 25, 30 25)), POLYGON ((25 28, 25 25, 20 25, 20 30, 25 28)), POLYGON ((30 30, 30 25, 25 25, 25 28, 30 30))) \ No newline at end of file