Skip to content

Add ST_Subdivide #719

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 1 commit into
base: main
Choose a base branch
from
Open

Add ST_Subdivide #719

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
87 changes: 87 additions & 0 deletions src/spatial/modules/geos/geos_geometry.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -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;
Expand Down Expand Up @@ -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;
Expand All @@ -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;
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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));
}
Expand Down Expand Up @@ -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);
}
Expand Down Expand Up @@ -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<size_t>(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<size_t>(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<size_t>(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 {
Expand Down
148 changes: 148 additions & 0 deletions src/spatial/modules/geos/geos_module.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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<double>::max();
xmax += std::numeric_limits<double>::max();
}

if (height == 0.0) {
ymin -= std::numeric_limits<double>::max();
ymax += std::numeric_limits<double>::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<string_t, uint32_t, string_t>(
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<ST_Touches> {
static bool ExecutePredicateNormal(const GeosGeometry &lhs, const GeosGeometry &rhs) {
return lhs.touches(rhs);
Expand Down Expand Up @@ -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);
Expand Down
41 changes: 41 additions & 0 deletions test/sql/geos/st_subdivide.test
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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)))