/********************************************************************** * * PostGIS - Spatial Types for PostgreSQL * http://postgis.net * * PostGIS is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 2 of the License, or * (at your option) any later version. * * PostGIS is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with PostGIS. If not, see . * ********************************************************************** * * Copyright 2019 Raúl Marín * **********************************************************************/ #include "lwgeom_wagyu.h" #define USE_WAGYU_INTERRUPT #include "mapbox/geometry/wagyu/wagyu.hpp" #include #include #include using namespace mapbox::geometry; using wagyu_coord_type = std::int32_t; using wagyu_polygon = mapbox::geometry::polygon; using wagyu_multipolygon = mapbox::geometry::multi_polygon; using wagyu_linearring = mapbox::geometry::linear_ring; using vwpolygon = std::vector; using wagyu_point = mapbox::geometry::point; static wagyu_linearring ptarray_to_wglinearring(const POINTARRAY *pa) { wagyu_linearring lr; lr.reserve(pa->npoints); size_t point_size = ptarray_point_size(pa); size_t pa_size = pa->npoints; uint8_t *buffer = pa->serialized_pointlist; for (std::uint32_t i = 0; i < pa_size; i++) { const wagyu_coord_type x = static_cast(*(double*) buffer); const wagyu_coord_type y = static_cast(((double*) buffer)[1]); buffer += point_size; lr.emplace_back(static_cast(x), static_cast(y)); } return lr; } static vwpolygon lwpoly_to_vwgpoly(const LWPOLY *geom, const GBOX *box) { LWCOLLECTION *geom_clipped_x = lwgeom_clip_to_ordinate_range((LWGEOM *)geom, 'X', box->xmin, box->xmax, 0); LWCOLLECTION *geom_clipped = lwgeom_clip_to_ordinate_range((LWGEOM *)geom_clipped_x, 'Y', box->ymin, box->ymax, 0); vwpolygon vp; for (std::uint32_t i = 0; i < geom_clipped->ngeoms; i++) { assert(geom_clipped->geoms[i]->type == POLYGONTYPE); LWPOLY *poly = (LWPOLY *)geom_clipped->geoms[i]; for (std::uint32_t ring = 0; ring < poly->nrings; ring++) { wagyu_polygon pol; pol.push_back(ptarray_to_wglinearring(poly->rings[ring])); /* If it has an inner ring, push it too */ ring++; if (ring != poly->nrings) { pol.push_back(ptarray_to_wglinearring(poly->rings[ring])); } vp.push_back(pol); } } lwgeom_free((LWGEOM *)geom_clipped_x); lwgeom_free((LWGEOM *)geom_clipped); return vp; } static vwpolygon lwmpoly_to_vwgpoly(const LWMPOLY *geom, const GBOX *box) { vwpolygon vp; for (std::uint32_t i = 0; i < geom->ngeoms; i++) { vwpolygon vp_intermediate = lwpoly_to_vwgpoly(geom->geoms[i], box); vp.insert(vp.end(), std::make_move_iterator(vp_intermediate.begin()), std::make_move_iterator(vp_intermediate.end())); } return vp; } /** * Transforms a liblwgeom geometry (types POLYGONTYPE or MULTIPOLYGONTYPE) * into an array of mapbox::geometry (polygon) * A single lwgeom polygon can lead to N mapbox polygon as odd numbered rings * are treated as new polygons (and even numbered treated as holes in the * previous rings) */ static vwpolygon lwgeom_to_vwgpoly(const LWGEOM *geom, const GBOX *box) { switch (geom->type) { case POLYGONTYPE: return lwpoly_to_vwgpoly(reinterpret_cast(geom), box); case MULTIPOLYGONTYPE: return lwmpoly_to_vwgpoly(reinterpret_cast(geom), box); default: return vwpolygon(); } } static POINTARRAY * wglinearring_to_ptarray(const wagyu_polygon::linear_ring_type &lr) { const uint32_t npoints = lr.size(); POINTARRAY *pa = ptarray_construct(0, 0, npoints); for (uint32_t i = 0; i < npoints; i++) { const wagyu_polygon::linear_ring_type::point_type &p = lr[i]; POINT4D point = {static_cast(p.x), static_cast(p.y), 0.0, 0.0}; ptarray_set_point4d(pa, i, &point); } return pa; } static LWGEOM * wgpoly_to_lwgeom(const wagyu_multipolygon::polygon_type &p) { const uint32_t nrings = p.size(); POINTARRAY **ppa = reinterpret_cast(lwalloc(sizeof(POINTARRAY *) * nrings)); for (uint32_t i = 0; i < nrings; i++) { ppa[i] = wglinearring_to_ptarray(p[i]); } return reinterpret_cast(lwpoly_construct(0, NULL, nrings, ppa)); } static LWGEOM * wgmpoly_to_lwgeom(const wagyu_multipolygon &mp) { const uint32_t ngeoms = mp.size(); if (ngeoms == 0) { return NULL; } else if (ngeoms == 1) { return wgpoly_to_lwgeom(mp[0]); } else { LWGEOM **geoms = reinterpret_cast(lwalloc(sizeof(LWGEOM *) * ngeoms)); for (uint32_t i = 0; i < ngeoms; i++) { geoms[i] = wgpoly_to_lwgeom(mp[i]); } return reinterpret_cast(lwcollection_construct(MULTIPOLYGONTYPE, 0, NULL, ngeoms, geoms)); } } static LWGEOM * __lwgeom_wagyu_clip_by_box(const LWGEOM *geom, const GBOX *bbox) { if (!geom || !bbox) { return NULL; } if (geom->type != POLYGONTYPE && geom->type != MULTIPOLYGONTYPE) { return NULL; } if (lwgeom_is_empty(geom)) { LWGEOM *out = lwgeom_construct_empty(MULTIPOLYGONTYPE, geom->srid, 0, 0); out->flags = geom->flags; return out; } wagyu_multipolygon solution; vwpolygon vpsubject = lwgeom_to_vwgpoly(geom, bbox); if (vpsubject.size() == 0) { LWGEOM *out = lwgeom_construct_empty(MULTIPOLYGONTYPE, geom->srid, 0, 0); out->flags = geom->flags; return out; } wagyu::wagyu clipper; for (auto &poly : vpsubject) { for (auto &lr : poly) { if (!lr.empty()) { clipper.add_ring(lr, wagyu::polygon_type_subject); } } } clipper.execute(wagyu::clip_type_union, solution, wagyu::fill_type_even_odd, wagyu::fill_type_even_odd); LWGEOM *g = wgmpoly_to_lwgeom(solution); if (g) g->srid = geom->srid; return g; } LWGEOM * lwgeom_wagyu_clip_by_box(const LWGEOM *geom, const GBOX *bbox) { mapbox::geometry::wagyu::interrupt_reset(); try { return __lwgeom_wagyu_clip_by_box(geom, bbox); } catch (...) { return NULL; } } const char * libwagyu_version() { static char str[50] = {0}; snprintf( str, sizeof(str), "%d.%d.%d (Internal)", WAGYU_MAJOR_VERSION, WAGYU_MINOR_VERSION, WAGYU_PATCH_VERSION); return str; } void lwgeom_wagyu_interruptRequest() { mapbox::geometry::wagyu::interrupt_request(); } void lwgeom_wagyu_interruptReset() { mapbox::geometry::wagyu::interrupt_reset(); }