Skip to content

Commit 6bd3397

Browse files
author
peng.li24
committed
feat: initial commit - C++ shapely library
Pixel-level C++ alignment of Python shapely geometry/ops functions. - geometry: Polygon, LineString, Point - ops: nearest_points, distance_to_multigeom
0 parents  commit 6bd3397

10 files changed

Lines changed: 871 additions & 0 deletions

‎geometry/linestring.h‎

Lines changed: 70 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,70 @@
1+
// Python Source: shapely/geometry/linestring.py
2+
// Line Range: class LineString(BaseGeometry)
3+
// Alignment: strict
4+
// EXEMPTION: cpp_template_optimization
5+
// Reason: C++ template to support both float32 and double coordinates.
6+
// Python shapely only has float64; float32 is a C++ zero-copy optimization.
7+
8+
#pragma once
9+
10+
#include <pybind11/pybind11.h>
11+
#include <pybind11/numpy.h>
12+
#include <memory>
13+
#include <geos/geom/GeometryFactory.h>
14+
#include <geos/geom/LineString.h>
15+
16+
namespace py = pybind11;
17+
18+
namespace shapely {
19+
namespace geometry {
20+
21+
#ifndef SHAPELY_GEOMETRY_POLYGON_DEFINED
22+
template <typename T>
23+
class Polygon;
24+
#endif
25+
26+
#ifndef SHAPELY_GEOMETRY_POINT_DEFINED
27+
template <typename T>
28+
class Point;
29+
#endif
30+
31+
template <typename T = double>
32+
class LineString {
33+
public:
34+
explicit LineString(const py::array_t<T>& coords);
35+
36+
LineString(LineString&&) = default;
37+
LineString& operator=(LineString&&) = default;
38+
39+
double distance(const LineString& other) const;
40+
template <typename U>
41+
double distance(const Polygon<U>& other) const;
42+
43+
template <typename U>
44+
double distance(const Point<U>& other) const;
45+
46+
template <typename U>
47+
bool intersects(const Polygon<U>& other) const;
48+
49+
template <typename U>
50+
double project(const Point<U>& other) const;
51+
52+
Point<double> interpolate(double distance) const;
53+
54+
double length() const;
55+
56+
py::array_t<T> coords_;
57+
58+
private:
59+
template <typename U> friend class Polygon;
60+
template <typename U> friend class Point;
61+
std::unique_ptr<geos::geom::LineString> geos_linestring_;
62+
geos::geom::GeometryFactory::Ptr factory_;
63+
};
64+
65+
} // namespace geometry
66+
} // namespace shapely
67+
68+
#define SHAPELY_GEOMETRY_LINESTRING_DEFINED
69+
70+
#include "linestring_impl.h"

‎geometry/linestring_impl.h‎

Lines changed: 102 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,102 @@
1+
// Python Source: shapely/geometry/linestring.py
2+
// Line Range: class LineString method implementations
3+
// Alignment: strict
4+
5+
#pragma once
6+
7+
#include "shapely/geometry/point.h"
8+
#include "shapely/geometry/polygon.h"
9+
10+
#include <geos/geom/GeometryFactory.h>
11+
#include <geos/geom/LineString.h>
12+
#include <geos/geom/Point.h>
13+
#include <geos/geom/Polygon.h>
14+
#include <geos/geom/Coordinate.h>
15+
#include <geos/geom/CoordinateSequence.h>
16+
#include <geos/geom/CoordinateSequenceFactory.h>
17+
#include <geos/operation/distance/DistanceOp.h>
18+
#include <geos/linearref/LengthIndexedLine.h>
19+
#include <stdexcept>
20+
21+
namespace shapely {
22+
namespace geometry {
23+
24+
template <typename T>
25+
LineString<T>::LineString(const py::array_t<T>& coords)
26+
: coords_(coords)
27+
{
28+
auto buf = coords.request();
29+
if (buf.ndim != 2 || buf.shape[1] < 2)
30+
throw std::runtime_error("LineString: coords must be [N, 2] or [N, 3] array");
31+
if (buf.shape[0] < 2)
32+
throw std::runtime_error("LineString: must have at least 2 points");
33+
34+
factory_ = geos::geom::GeometryFactory::create();
35+
T *p = static_cast<T *>(buf.ptr);
36+
size_t n = buf.shape[0];
37+
int stride = buf.shape[1];
38+
39+
auto coord_seq = factory_->getCoordinateSequenceFactory()->create(n, 2);
40+
for (size_t i = 0; i < n; ++i)
41+
coord_seq->setAt(geos::geom::Coordinate(
42+
static_cast<double>(p[i * stride + 0]),
43+
static_cast<double>(p[i * stride + 1])), i);
44+
45+
geos_linestring_ = factory_->createLineString(std::move(coord_seq));
46+
}
47+
48+
template <typename T>
49+
double LineString<T>::distance(const LineString& other) const
50+
{
51+
geos::operation::distance::DistanceOp dist_op(geos_linestring_.get(), other.geos_linestring_.get());
52+
return dist_op.distance();
53+
}
54+
55+
template <typename T>
56+
template <typename U>
57+
double LineString<T>::distance(const Polygon<U>& other) const
58+
{
59+
geos::operation::distance::DistanceOp dist_op(geos_linestring_.get(), other.geos_polygon_.get());
60+
return dist_op.distance();
61+
}
62+
63+
template <typename T>
64+
template <typename U>
65+
double LineString<T>::distance(const Point<U>& other) const
66+
{
67+
geos::operation::distance::DistanceOp dist_op(geos_linestring_.get(), other.geos_point_.get());
68+
return dist_op.distance();
69+
}
70+
71+
template <typename T>
72+
template <typename U>
73+
bool LineString<T>::intersects(const Polygon<U>& other) const
74+
{
75+
return geos_linestring_->intersects(other.geos_polygon_.get());
76+
}
77+
78+
template <typename T>
79+
template <typename U>
80+
double LineString<T>::project(const Point<U>& other) const
81+
{
82+
geos::linearref::LengthIndexedLine indexed_line(geos_linestring_.get());
83+
return indexed_line.project(geos::geom::Coordinate(
84+
static_cast<double>(other.x_), static_cast<double>(other.y_)));
85+
}
86+
87+
template <typename T>
88+
Point<double> LineString<T>::interpolate(double distance) const
89+
{
90+
geos::linearref::LengthIndexedLine indexed_line(geos_linestring_.get());
91+
auto coord = indexed_line.extractPoint(distance);
92+
return Point<double>(coord.x, coord.y);
93+
}
94+
95+
template <typename T>
96+
double LineString<T>::length() const
97+
{
98+
return geos_linestring_->getLength();
99+
}
100+
101+
} // namespace geometry
102+
} // namespace shapely

‎geometry/point.h‎

Lines changed: 63 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,63 @@
1+
// Python Source: shapely/geometry/point.py
2+
// Line Range: class Point(BaseGeometry)
3+
// Alignment: strict
4+
// EXEMPTION: cpp_template_optimization
5+
// Reason: C++ template to support both float32 and double coordinates.
6+
// Python shapely only has float64; float32 is a C++ zero-copy optimization.
7+
8+
#pragma once
9+
10+
#include <pybind11/pybind11.h>
11+
#include <pybind11/numpy.h>
12+
#include <memory>
13+
#include <geos/geom/GeometryFactory.h>
14+
#include <geos/geom/Point.h>
15+
16+
namespace py = pybind11;
17+
18+
namespace shapely {
19+
namespace geometry {
20+
21+
#ifndef SHAPELY_GEOMETRY_LINESTRING_DEFINED
22+
template <typename T>
23+
class LineString;
24+
#endif
25+
26+
#ifndef SHAPELY_GEOMETRY_POLYGON_DEFINED
27+
template <typename T>
28+
class Polygon;
29+
#endif
30+
31+
template <typename T = double>
32+
class Point {
33+
public:
34+
Point(T x, T y);
35+
36+
Point(Point&&) = default;
37+
Point& operator=(Point&&) = default;
38+
39+
double distance(const Point& other) const;
40+
41+
template <typename U>
42+
double distance(const LineString<U>& other) const;
43+
44+
template <typename U>
45+
double distance(const Polygon<U>& other) const;
46+
47+
Polygon<double> buffer(double distance) const;
48+
49+
T x_, y_;
50+
51+
private:
52+
template <typename U> friend class LineString;
53+
template <typename U> friend class Polygon;
54+
std::unique_ptr<geos::geom::Point> geos_point_;
55+
geos::geom::GeometryFactory::Ptr factory_;
56+
};
57+
58+
} // namespace geometry
59+
} // namespace shapely
60+
61+
#define SHAPELY_GEOMETRY_POINT_DEFINED
62+
63+
#include "point_impl.h"

‎geometry/point_impl.h‎

Lines changed: 92 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,92 @@
1+
// Python Source: shapely/geometry/point.py
2+
// Line Range: class Point method implementations
3+
// Alignment: strict
4+
5+
#pragma once
6+
7+
#include "shapely/geometry/linestring.h"
8+
#include "shapely/geometry/polygon.h"
9+
10+
#include <geos/geom/GeometryFactory.h>
11+
#include <geos/geom/Point.h>
12+
#include <geos/geom/LineString.h>
13+
#include <geos/geom/Polygon.h>
14+
#include <geos/geom/Coordinate.h>
15+
#include <geos/geom/CoordinateSequence.h>
16+
#include <geos/geom/CoordinateSequenceFactory.h>
17+
#include <geos/operation/distance/DistanceOp.h>
18+
#include <stdexcept>
19+
20+
namespace shapely {
21+
namespace geometry {
22+
23+
template <typename T>
24+
Point<T>::Point(T x, T y)
25+
: x_(x), y_(y)
26+
{
27+
factory_ = geos::geom::GeometryFactory::create();
28+
geos_point_.reset(factory_->createPoint(
29+
geos::geom::Coordinate(static_cast<double>(x), static_cast<double>(y))));
30+
}
31+
32+
template <typename T>
33+
double Point<T>::distance(const Point& other) const
34+
{
35+
geos::operation::distance::DistanceOp dist_op(geos_point_.get(), other.geos_point_.get());
36+
return dist_op.distance();
37+
}
38+
39+
template <typename T>
40+
template <typename U>
41+
double Point<T>::distance(const LineString<U>& other) const
42+
{
43+
geos::operation::distance::DistanceOp dist_op(geos_point_.get(), other.geos_linestring_.get());
44+
return dist_op.distance();
45+
}
46+
47+
template <typename T>
48+
template <typename U>
49+
double Point<T>::distance(const Polygon<U>& other) const
50+
{
51+
geos::operation::distance::DistanceOp dist_op(geos_point_.get(), other.geos_polygon_.get());
52+
return dist_op.distance();
53+
}
54+
55+
template <typename T>
56+
Polygon<double> Point<T>::buffer(double distance) const
57+
{
58+
auto buf_geom = geos_point_->buffer(distance);
59+
if (buf_geom == nullptr || buf_geom->isEmpty())
60+
return std::move(Polygon<double>(py::array_t<double>(std::vector<py::ssize_t>{0, 2})));
61+
62+
const geos::geom::Geometry *poly = buf_geom.get();
63+
if (poly->getGeometryTypeId() != geos::geom::GEOS_POLYGON)
64+
{
65+
if (poly->getNumGeometries() > 0)
66+
poly = poly->getGeometryN(0);
67+
}
68+
69+
if (poly->getGeometryTypeId() != geos::geom::GEOS_POLYGON || poly->isEmpty())
70+
return std::move(Polygon<double>(py::array_t<double>(std::vector<py::ssize_t>{0, 2})));
71+
72+
const geos::geom::Polygon *geos_poly = dynamic_cast<const geos::geom::Polygon *>(poly);
73+
if (!geos_poly)
74+
return std::move(Polygon<double>(py::array_t<double>(std::vector<py::ssize_t>{0, 2})));
75+
76+
const geos::geom::CoordinateSequence *cs = geos_poly->getExteriorRing()->getCoordinatesRO();
77+
if (!cs || cs->isEmpty())
78+
return std::move(Polygon<double>(py::array_t<double>(std::vector<py::ssize_t>{0, 2})));
79+
80+
size_t n = cs->getSize();
81+
py::array_t<double> coords(std::vector<py::ssize_t>{static_cast<py::ssize_t>(n - 1), static_cast<py::ssize_t>(2)});
82+
double *p = static_cast<double *>(coords.request().ptr);
83+
for (size_t i = 0; i < n - 1; ++i) {
84+
p[i * 2] = cs->getAt(i).x;
85+
p[i * 2 + 1] = cs->getAt(i).y;
86+
}
87+
88+
return std::move(Polygon<double>(coords));
89+
}
90+
91+
} // namespace geometry
92+
} // namespace shapely

‎geometry/polygon.h‎

Lines changed: 67 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,67 @@
1+
// Python Source: shapely/geometry/polygon.py
2+
// Line Range: class Polygon(BaseGeometry)
3+
// Alignment: strict
4+
// EXEMPTION: cpp_template_optimization
5+
// Reason: C++ template to support both float32 and double coordinates.
6+
// Python shapely only has float64; float32 is a C++ zero-copy optimization.
7+
8+
#pragma once
9+
10+
#include <pybind11/pybind11.h>
11+
#include <pybind11/numpy.h>
12+
#include <memory>
13+
#include <geos/geom/GeometryFactory.h>
14+
#include <geos/geom/Polygon.h>
15+
16+
namespace py = pybind11;
17+
18+
namespace shapely {
19+
namespace geometry {
20+
21+
#ifndef SHAPELY_GEOMETRY_LINESTRING_DEFINED
22+
template <typename T>
23+
class LineString;
24+
#endif
25+
26+
template <typename T = double>
27+
class Polygon {
28+
public:
29+
explicit Polygon(const py::array_t<T>& coords);
30+
31+
Polygon(Polygon&&) = default;
32+
Polygon& operator=(Polygon&&) = default;
33+
34+
double area() const;
35+
double distance(const Polygon& other) const;
36+
37+
template <typename U>
38+
double distance(const LineString<U>& other) const;
39+
40+
bool intersects(const Polygon& other) const;
41+
42+
template <typename U>
43+
bool intersects(const LineString<U>& other) const;
44+
45+
Polygon<double> intersection(const Polygon<double>& other) const;
46+
double intersection_area(const Polygon<double>& other) const;
47+
bool is_valid() const;
48+
bool is_empty() const;
49+
Polygon<double> buffer(double distance) const;
50+
py::array_t<T> exterior_coords() const;
51+
52+
py::array_t<T> coords_;
53+
54+
private:
55+
template <typename U> friend class Polygon;
56+
template <typename U> friend class LineString;
57+
template <typename U> friend class Point;
58+
std::unique_ptr<geos::geom::Polygon> geos_polygon_;
59+
geos::geom::GeometryFactory::Ptr factory_;
60+
};
61+
62+
} // namespace geometry
63+
} // namespace shapely
64+
65+
#define SHAPELY_GEOMETRY_POLYGON_DEFINED
66+
67+
#include "polygon_impl.h"

0 commit comments

Comments
 (0)