Skip to content

Commit

Permalink
Predicates functor (#22)
Browse files Browse the repository at this point in the history 
* clang-format

* add Predicates functor

A bit more DRY. For now also allows creating the
S2BooleanOperator::Options object only once (although we might want to
eventually expose options as optional arguments in Python).

* fix compilation warnings

Unrelated to this PR.

* use functor for all implemented predicates
  • Loading branch information
@benbovy
benbovy authored Mar 13, 2023
1 parent e8bf2e5 commit eb437b6
Show file tree
Hide file tree
Showing 4 changed files with 89 additions and 84 deletions.
63 changes: 30 additions & 33 deletions src/geography.cpp
View file
Original file line number Diff line number Diff line change
@@ -1,13 +1,5 @@
#include "geography.hpp"

#include <memory>
#include <stdexcept>
#include <sstream>
#include <type_traits>
#include <utility>
#include <variant>
#include <vector>

#include <pybind11/pybind11.h>
#include <pybind11/stl.h>
#include <s2/s2latlng.h>
Expand All @@ -16,6 +8,14 @@
#include <s2/s2polygon.h>
#include <s2geography.h>

#include <memory>
#include <sstream>
#include <stdexcept>
#include <type_traits>
#include <utility>
#include <variant>
#include <vector>

#include "pybind11.hpp"
#include "pybind11/stl.h"

Expand All @@ -38,13 +38,11 @@ py::detail::type_info *PyObjectGeography::geography_tinfo = nullptr;
// It is not fully supported with Pybind11 (Point is non-default constructible)
// https://github.com/pybind/pybind11/issues/4108
//
S2Point to_s2point(const std::pair<double, double>& vertex) {
S2Point to_s2point(const std::pair<double, double> &vertex) {
return S2LatLng::FromDegrees(vertex.first, vertex.second).ToPoint();
}

S2Point to_s2point(const Point* vertex) {
return vertex->s2point();
}
S2Point to_s2point(const Point *vertex) { return vertex->s2point(); }

/*
** Helper to create Geography object wrappers.
Expand All @@ -54,12 +52,11 @@ S2Point to_s2point(const Point* vertex) {
** @tparam S The S2Geometry type
*/
template <class T1, class T2, class S>
std::unique_ptr<T2> make_geography(S&& s2_obj) {
std::unique_ptr<T2> make_geography(S &&s2_obj) {
S2GeographyPtr s2geog_ptr = std::make_unique<T1>(std::forward<S>(s2_obj));
return std::make_unique<T2>(std::move(s2geog_ptr));
}


class PointFactory {
public:
static std::unique_ptr<Point> FromLatLonDegrees(double lat_degrees,
Expand All @@ -69,29 +66,30 @@ class PointFactory {
return make_geography<s2geog::PointGeography, Point>(S2Point(latlng));
}

//TODO: from LatLonRadians
// TODO: from LatLonRadians
};

template <class V>
static std::unique_ptr<LineString> create_linestring(const std::vector<V>& coords) {
static std::unique_ptr<LineString> create_linestring(
const std::vector<V> &coords) {
std::vector<S2Point> pts(coords.size());

std::transform(
coords.begin(), coords.end(), pts.begin(),
[](const V& vertex) { return to_s2point(vertex); });
std::transform(coords.begin(), coords.end(), pts.begin(),
[](const V &vertex) { return to_s2point(vertex); });

auto polyline_ptr = std::make_unique<S2Polyline>(pts);

return make_geography<s2geog::PolylineGeography, LineString>(std::move(polyline_ptr));
return make_geography<s2geog::PolylineGeography, LineString>(
std::move(polyline_ptr));
}

template <class V>
static std::unique_ptr<spherely::Polygon> create_polygon(const std::vector<V>& shell) {
static std::unique_ptr<spherely::Polygon> create_polygon(
const std::vector<V> &shell) {
std::vector<S2Point> shell_pts(shell.size());

std::transform(
shell.begin(), shell.end(), shell_pts.begin(),
[](const V& vertex) { return to_s2point(vertex); });
std::transform(shell.begin(), shell.end(), shell_pts.begin(),
[](const V &vertex) { return to_s2point(vertex); });

auto shell_loop_ptr = std::make_unique<S2Loop>();
// TODO: maybe add an option to skip validity checks
Expand All @@ -118,7 +116,8 @@ static std::unique_ptr<spherely::Polygon> create_polygon(const std::vector<V>& s
polygon_ptr->set_s2debug_override(S2Debug::DISABLE);
polygon_ptr->InitOriented(std::move(loops));

return make_geography<s2geog::PolygonGeography, spherely::Polygon>(std::move(polygon_ptr));
return make_geography<s2geog::PolygonGeography, spherely::Polygon>(
std::move(polygon_ptr));
}

/*
Expand All @@ -132,7 +131,7 @@ py::array_t<int> num_shapes(const py::array_t<PyObjectGeography> geographies) {
py::buffer_info result_buf = result.request();
int *rptr = static_cast<int *>(result_buf.ptr);

for (size_t i = 0; i < buf.size; i++) {
for (py::ssize_t i = 0; i < buf.size; i++) {
auto geog_ptr = (*geographies.data(i)).as_geog_ptr();
rptr[i] = geog_ptr->num_shapes();
}
Expand All @@ -157,9 +156,7 @@ py::array_t<PyObjectGeography> create(py::array_t<double> xs,
double *yptr = static_cast<double *>(ybuf.ptr);
py::object *rptr = static_cast<py::object *>(rbuf.ptr);

size_t size = static_cast<size_t>(xbuf.shape[0]);

for (size_t i = 0; i < xbuf.shape[0]; i++) {
for (py::ssize_t i = 0; i < xbuf.shape[0]; i++) {
auto point_ptr = PointFactory::FromLatLonDegrees(xptr[i], yptr[i]);
// rptr[i] = PyObjectGeography::as_py_object(std::move(point_ptr));
rptr[i] = py::cast(std::move(point_ptr));
Expand Down Expand Up @@ -283,10 +280,9 @@ void init_geography(py::module &m) {
pylinestring.def(py::init(&create_linestring<std::pair<double, double>>),
py::arg("coordinates"));

pylinestring.def(py::init(&create_linestring<Point*>),
pylinestring.def(py::init(&create_linestring<Point *>),
py::arg("coordinates"));


auto pypolygon =
py::class_<spherely::Polygon, Geography>(m, "Polygon", R"pbdoc(
A geography type representing an area that is enclosed by a linear ring.
Expand All @@ -301,8 +297,9 @@ void init_geography(py::module &m) {
)pbdoc");

pypolygon.def(py::init(&create_polygon<std::pair<double, double>>), py::arg("shell"));
pypolygon.def(py::init(&create_polygon<Point*>), py::arg("shell"));
pypolygon.def(py::init(&create_polygon<std::pair<double, double>>),
py::arg("shell"));
pypolygon.def(py::init(&create_polygon<Point *>), py::arg("shell"));

// Temp test

Expand Down
10 changes: 8 additions & 2 deletions src/geography.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,12 @@ using S2GeographyIndexPtr = std::unique_ptr<s2geog::ShapeIndexGeography>;
/*
** The registered Geography types
*/
enum class GeographyType : std::int8_t { None = -1, Point, LineString, Polygon };
enum class GeographyType : std::int8_t {
None = -1,
Point,
LineString,
Polygon
};

/*
** Thin wrapper around s2geography::Geography.
Expand Down Expand Up @@ -81,7 +86,8 @@ class Point : public Geography {
}

inline const S2Point& s2point() const {
const auto& points = static_cast<const s2geog::PointGeography&>(geog()).Points();
const auto& points =
static_cast<const s2geog::PointGeography&>(geog()).Points();
// TODO: does not work for empty point geography
return points[0];
}
Expand Down
97 changes: 49 additions & 48 deletions src/predicates.cpp
View file
Original file line number Diff line number Diff line change
@@ -1,55 +1,41 @@
#include <s2/s2boolean_operation.h>
#include <s2geography.h>

#include <functional>

#include "geography.hpp"
#include "pybind11.hpp"

namespace py = pybind11;
namespace s2geog = s2geography;
using namespace spherely;

bool intersects(PyObjectGeography a, PyObjectGeography b) {
const auto& a_index = a.as_geog_ptr()->geog_index();
const auto& b_index = b.as_geog_ptr()->geog_index();

S2BooleanOperation::Options options;
return s2geog::s2_intersects(a_index, b_index, options);
}

bool equals(PyObjectGeography a, PyObjectGeography b) {
const auto& a_index = a.as_geog_ptr()->geog_index();
const auto& b_index = b.as_geog_ptr()->geog_index();

S2BooleanOperation::Options options;
return s2geog::s2_equals(a_index, b_index, options);
}

bool contains(PyObjectGeography a, PyObjectGeography b) {
const auto& a_index = a.as_geog_ptr()->geog_index();
const auto& b_index = b.as_geog_ptr()->geog_index();

S2BooleanOperation::Options options;
return s2geog::s2_contains(a_index, b_index, options);
}

bool within(PyObjectGeography a, PyObjectGeography b) {
const auto& a_index = a.as_geog_ptr()->geog_index();
const auto& b_index = b.as_geog_ptr()->geog_index();

S2BooleanOperation::Options options;
return s2geog::s2_contains(b_index, a_index, options);
}
/*
** Functor for predicate bindings.
*/
class Predicate {
public:
using FuncType = std::function<bool(const s2geog::ShapeIndexGeography&,
const s2geog::ShapeIndexGeography&,
const S2BooleanOperation::Options&)>;

template <class F>
Predicate(F&& func) : m_func(std::forward<F>(func)) {}

bool operator()(PyObjectGeography a, PyObjectGeography b) const {
const auto& a_index = a.as_geog_ptr()->geog_index();
const auto& b_index = b.as_geog_ptr()->geog_index();
return m_func(a_index, b_index, m_options);
}

bool disjoint(PyObjectGeography a, PyObjectGeography b) {
const auto& a_index = a.as_geog_ptr()->geog_index();
const auto& b_index = b.as_geog_ptr()->geog_index();

S2BooleanOperation::Options options;
return !s2geog::s2_intersects(a_index, b_index, options);
}
private:
FuncType m_func;
S2BooleanOperation::Options m_options;
};

void init_predicates(py::module& m) {
m.def("intersects", py::vectorize(&intersects), py::arg("a"), py::arg("b"),
m.def("intersects", py::vectorize(Predicate(s2geog::s2_intersects)),
py::arg("a"), py::arg("b"),
R"pbdoc(
Returns True if A and B share any portion of space.
Expand All @@ -62,7 +48,8 @@ void init_predicates(py::module& m) {
)pbdoc");

m.def("equals", py::vectorize(&equals), py::arg("a"), py::arg("b"),
m.def("equals", py::vectorize(Predicate(s2geog::s2_equals)), py::arg("a"),
py::arg("b"),
R"pbdoc(
Returns True if A and B are spatially equal.
Expand All @@ -76,7 +63,8 @@ void init_predicates(py::module& m) {
)pbdoc");

m.def("contains", py::vectorize(&contains), py::arg("a"), py::arg("b"),
m.def("contains", py::vectorize(Predicate(s2geog::s2_contains)),
py::arg("a"), py::arg("b"),
R"pbdoc(
Returns True if B is completely inside A.
Expand All @@ -86,9 +74,16 @@ void init_predicates(py::module& m) {
Geography object(s)
)pbdoc");

m.def("within", py::vectorize(&within), py::arg("a"), py::arg("b"),
R"pbdoc(

m.def(
"within",
py::vectorize(Predicate([](const s2geog::ShapeIndexGeography& a_index,
const s2geog::ShapeIndexGeography& b_index,
const S2BooleanOperation::Options& options) {
return s2geog::s2_contains(b_index, a_index, options);
})),
py::arg("a"), py::arg("b"),
R"pbdoc(
Returns True if A is completely inside B.
Parameters
Expand All @@ -98,8 +93,15 @@ void init_predicates(py::module& m) {
)pbdoc");

m.def("disjoint", py::vectorize(&disjoint), py::arg("a"), py::arg("b"),
R"pbdoc(
m.def(
"disjoint",
py::vectorize(Predicate([](const s2geog::ShapeIndexGeography& a_index,
const s2geog::ShapeIndexGeography& b_index,
const S2BooleanOperation::Options& options) {
return !s2geog::s2_intersects(a_index, b_index, options);
})),
py::arg("a"), py::arg("b"),
R"pbdoc(
Returns True if A boundaries and interior does not intersect at all
with those of B.
Expand All @@ -109,5 +111,4 @@ void init_predicates(py::module& m) {
Geography object(s)
)pbdoc");

}
}
3 changes: 2 additions & 1 deletion src/pybind11.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -142,7 +142,8 @@ struct npy_format_descriptor<spherely::PyObjectGeography> {
};

// Override signature type hint for vectorized Geography arguments
template <int Flags> struct handle_type_name<array_t<spherely::PyObjectGeography, Flags>> {
template <int Flags>
struct handle_type_name<array_t<spherely::PyObjectGeography, Flags>> {
static constexpr auto name = _("Geography | array_like");
};

Expand Down

0 comments on commit eb437b6

Please sign in to comment.