Skip to content

Commit 9ded4e5

Browse files
author
peng.li24
committed
feat: update geometry headers, module, and tests
1 parent 1dadc54 commit 9ded4e5

8 files changed

Lines changed: 118 additions & 140 deletions

File tree

‎pycpp/geometry_py.h‎

Lines changed: 5 additions & 48 deletions
Original file line numberDiff line numberDiff line change
@@ -22,7 +22,6 @@
2222
#include "../shapely/geometry/geometrycollection.h"
2323

2424
#include <array>
25-
#include <cmath>
2625
#include <cstring>
2726
#include <tuple>
2827
#include <vector>
@@ -67,20 +66,6 @@ inline void ensure_c_contiguous(const py::buffer_info& buf) {
6766
}
6867
}
6968

70-
/// Validate coordinate buffer — reject NaN/Inf before they reach GEOS.
71-
/// Throws std::invalid_argument with the index of the first offending value.
72-
template <typename T>
73-
inline void ensure_no_nan_inf(const T* data, size_t count, const char* label = "coords") {
74-
for (size_t i = 0; i < count; ++i) {
75-
if (std::isnan(data[i]))
76-
throw std::invalid_argument(
77-
std::string(label) + " contains NaN at index " + std::to_string(i));
78-
if (std::isinf(data[i]))
79-
throw std::invalid_argument(
80-
std::string(label) + " contains Inf at index " + std::to_string(i));
81-
}
82-
}
83-
8469
/// Convert any py::array coordinate buffer to std::vector<double>.
8570
/// Extracts the first 2 columns per row — safe for (N,2) and (N,≥3) arrays.
8671
/// dtype double → direct copy; float32/other → cast & widen to double.
@@ -111,8 +96,6 @@ inline std::vector<double> array_to_double_vec(const py::array& arr) {
11196
tmp[r * 2 + 1] = static_cast<double>(src[base + 1]);
11297
}
11398
}
114-
// Validate AFTER conversion: reject NaN/Inf before GEOS
115-
ensure_no_nan_inf(tmp.data(), rows * 2, "coords");
11699
return tmp;
117100
}
118101

@@ -121,65 +104,45 @@ inline std::vector<double> array_to_double_vec(const py::array& arr) {
121104
// ============================================================================
122105

123106
// -- Point: scalar (x,y) + array ([x,y]) forms --
124-
inline Point<double> point(double x, double y) {
125-
if (std::isnan(x) || std::isnan(y))
126-
throw std::invalid_argument("point coords contains NaN");
127-
if (std::isinf(x) || std::isinf(y))
128-
throw std::invalid_argument("point coords contains Inf");
129-
return Point<double>(x, y);
130-
}
107+
inline Point<double> point(double x, double y) { return Point<double>(x, y); }
131108
inline Point<float> point(float x, float y) { return Point<float>(x, y); }
132109

133110
inline Point<double> point(const py::array_t<double>& arr) {
134111
auto buf = arr.request();
135112
ensure_c_contiguous(buf);
136113
const double* p = static_cast<const double*>(buf.ptr);
137-
ensure_no_nan_inf(p, buf.size, "point coords");
138114
return Point<double>(buf.size > 0 ? p[0] : 0, buf.size > 1 ? p[1] : 0);
139115
}
140116
// float32 auto-upcast: prevent GEOS precision mismatch vs Python shapely
141117
inline Point<double> point(const py::array_t<float>& arr) {
142118
auto buf = arr.request();
143119
ensure_c_contiguous(buf);
144120
const float* p = static_cast<const float*>(buf.ptr);
145-
double x = buf.size > 0 ? static_cast<double>(p[0]) : 0;
146-
double y = buf.size > 1 ? static_cast<double>(p[1]) : 0;
147-
if (std::isnan(x) || std::isnan(y))
148-
throw std::invalid_argument("point coords contains NaN");
149-
if (std::isinf(x) || std::isinf(y))
150-
throw std::invalid_argument("point coords contains Inf");
151-
return Point<double>(x, y);
121+
return Point<double>(buf.size > 0 ? static_cast<double>(p[0]) : 0,
122+
buf.size > 1 ? static_cast<double>(p[1]) : 0);
152123
}
153124
// auto-double: accept any dtype (int, unknown), always produce Point<double>
154125
// Handles 1D [x,y] and 2D (N,≥2) arrays — takes first 2 elements.
155126
inline Point<double> point(const py::array& arr) {
156127
auto buf = arr.request();
157128
ensure_c_contiguous(buf);
158129
if (buf.size < 2) return Point<double>(0, 0);
159-
double x, y;
160130
if (arr.dtype().is(py::dtype::of<double>())) {
161131
const double* p = static_cast<const double*>(buf.ptr);
162-
x = p[0]; y = p[1];
132+
return Point<double>(p[0], p[1]);
163133
} else {
164134
// float32 or int/other → cast through pybind11 to float32, then widen
165135
auto f32 = py::cast<py::array_t<float>>(arr);
166136
auto fbuf = f32.request();
167137
const float* p = static_cast<const float*>(fbuf.ptr);
168-
x = static_cast<double>(p[0]);
169-
y = static_cast<double>(p[1]);
138+
return Point<double>(static_cast<double>(p[0]), static_cast<double>(p[1]));
170139
}
171-
if (std::isnan(x) || std::isnan(y))
172-
throw std::invalid_argument("point coords contains NaN");
173-
if (std::isinf(x) || std::isinf(y))
174-
throw std::invalid_argument("point coords contains Inf");
175-
return Point<double>(x, y);
176140
}
177141

178142
// -- LineString --
179143
inline LineString<double> linestring(const py::array_t<double>& arr) {
180144
auto buf = arr.request();
181145
ensure_c_contiguous(buf);
182-
ensure_no_nan_inf(static_cast<const double*>(buf.ptr), buf.size, "linestring coords");
183146
return LineString<double>(static_cast<const double*>(buf.ptr), buf.shape[0], buf.shape[1]);
184147
}
185148
// float32 auto-upcast: prevent GEOS precision mismatch vs Python shapely
@@ -204,15 +167,13 @@ inline LineString<double> linestring(const std::vector<std::array<double, 2>>& p
204167
flat[i * 2] = pts[i][0];
205168
flat[i * 2 + 1] = pts[i][1];
206169
}
207-
ensure_no_nan_inf(flat.data(), pts.size() * 2, "linestring coords");
208170
return LineString<double>(flat.data(), pts.size(), 2);
209171
}
210172

211173
// -- Polygon --
212174
inline Polygon<double> polygon(const py::array_t<double>& arr) {
213175
auto buf = arr.request();
214176
ensure_c_contiguous(buf);
215-
ensure_no_nan_inf(static_cast<const double*>(buf.ptr), buf.size, "polygon coords");
216177
return Polygon<double>(static_cast<const double*>(buf.ptr), buf.shape[0], buf.shape[1]);
217178
}
218179
// float32 auto-upcast: prevent GEOS precision mismatch vs Python shapely
@@ -231,7 +192,6 @@ inline Polygon<double> polygon(const py::array& arr) {
231192
inline LinearRing<double> linearring(const py::array_t<double>& arr) {
232193
auto buf = arr.request();
233194
ensure_c_contiguous(buf);
234-
ensure_no_nan_inf(static_cast<const double*>(buf.ptr), buf.size, "linearring coords");
235195
return LinearRing<double>(static_cast<const double*>(buf.ptr), buf.shape[0], buf.shape[1]);
236196
}
237197
inline LinearRing<double> linearring(const py::array& arr) {
@@ -244,7 +204,6 @@ inline LinearRing<double> linearring(const py::array& arr) {
244204
inline MultiPoint<double> multipoint(const py::array_t<double>& arr) {
245205
auto buf = arr.request();
246206
ensure_c_contiguous(buf);
247-
ensure_no_nan_inf(static_cast<const double*>(buf.ptr), buf.size, "multipoint coords");
248207
return MultiPoint<double>(static_cast<const double*>(buf.ptr), buf.shape[0], buf.shape[1]);
249208
}
250209
// float32 auto-upcast: prevent GEOS precision mismatch vs Python shapely
@@ -265,7 +224,6 @@ inline MultiLineString<double> multilinestring(const std::vector<py::array_t<dou
265224
for (auto& arr : arrays) {
266225
auto buf = arr.request();
267226
ensure_c_contiguous(buf);
268-
ensure_no_nan_inf(static_cast<const double*>(buf.ptr), buf.size, "multilinestring coords");
269227
mls.add_line(static_cast<const double*>(buf.ptr), buf.shape[0], buf.shape[1]);
270228
}
271229
return mls;
@@ -296,7 +254,6 @@ inline MultiPolygon<double> multipolygon(const std::vector<py::array_t<double>>&
296254
for (auto& arr : arrays) {
297255
auto buf = arr.request();
298256
ensure_c_contiguous(buf);
299-
ensure_no_nan_inf(static_cast<const double*>(buf.ptr), buf.size, "multipolygon coords");
300257
mp.add_polygon(static_cast<const double*>(buf.ptr), buf.shape[0], buf.shape[1]);
301258
}
302259
return mp;

‎shapely/geometry/base.h‎

Lines changed: 41 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -31,9 +31,34 @@
3131
#include <geos/simplify/TopologyPreservingSimplifier.h>
3232
#include <geos/operation/valid/IsValidOp.h>
3333

34+
#include <cmath>
35+
#include <stdexcept>
36+
3437
namespace shapely {
3538
namespace detail {
3639

40+
/// Catch GEOS exceptions and return a safe fallback (matching Python shapely behavior).
41+
/// Python shapely does NOT throw on NaN/Inf geometries — it returns NaN / inf / False.
42+
#define GEOS_SAFE_DOUBLE(op_name, expr) \
43+
([&]() -> double { \
44+
try { return (expr); } \
45+
catch (const std::exception&) { \
46+
return std::numeric_limits<double>::quiet_NaN(); \
47+
} \
48+
}())
49+
50+
#define GEOS_SAFE_BOOL(op_name, expr) \
51+
([&]() -> bool { \
52+
try { return (expr); } \
53+
catch (const std::exception&) { return false; } \
54+
}())
55+
56+
#define GEOS_SAFE_STRING(op_name, expr) \
57+
([&]() -> std::string { \
58+
try { return (expr); } \
59+
catch (const std::exception&) { return std::string(); } \
60+
}())
61+
3762
// --- WKT / WKB serialization ------------------------------------------------
3863

3964
// Python: shapely/geometry/base.py::wkt:L369
@@ -194,14 +219,27 @@ inline double geos_length(const geos::geom::Geometry* g) {
194219
// Python: shapely/geometry/base.py::distance:L438
195220
inline double geos_distance(const geos::geom::Geometry* a, const geos::geom::Geometry* b) {
196221
if (!a || !b) return 0.0;
197-
geos::operation::distance::DistanceOp op(a, b);
198-
return op.distance();
222+
return GEOS_SAFE_DOUBLE("distance", [&]{
223+
geos::operation::distance::DistanceOp op(a, b);
224+
return op.distance();
225+
}());
199226
}
200227

201228
// Python: shapely/geometry/base.py::hausdorff_distance:L442
202229
inline double geos_hausdorff_distance(const geos::geom::Geometry* a, const geos::geom::Geometry* b) {
203230
if (!a || !b) return 0.0;
204-
return geos::algorithm::distance::DiscreteHausdorffDistance(*a, *b).distance();
231+
return GEOS_SAFE_DOUBLE("hausdorff_distance",
232+
geos::algorithm::distance::DiscreteHausdorffDistance(*a, *b).distance());
233+
}
234+
235+
/// Safe DistanceOp wrapper — returns NaN on GEOS exception (matching Python shapely).
236+
/// Python shapely does NOT crash on degenerate/NaN geometries; it returns NaN/inf.
237+
inline double geos_safe_distance(const geos::geom::Geometry* a, const geos::geom::Geometry* b) {
238+
if (!a || !b) return 0.0;
239+
return GEOS_SAFE_DOUBLE("distance", [&] {
240+
geos::operation::distance::DistanceOp op(a, b);
241+
return op.distance();
242+
}());
205243
}
206244

207245
// Python: shapely/geometry/base.py::bounds:L470

‎shapely/geometry/linestring.h‎

Lines changed: 4 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -234,22 +234,19 @@ std::tuple<std::vector<T>, std::vector<T>> LineString<T>::xy() const {
234234

235235
template <typename T>
236236
double LineString<T>::distance(const LineString& o) const {
237-
geos::operation::distance::DistanceOp op(geos_linestring_.get(), o.geos_linestring_.get());
238-
return op.distance();
237+
return detail::geos_safe_distance(geos_linestring_.get(), o.geos_linestring_.get());
239238
}
240239
template <typename T> template <typename U>
241240
double LineString<T>::distance(const Polygon<U>& o) const {
242-
geos::operation::distance::DistanceOp op(geos_linestring_.get(), o.geos_polygon_.get());
243-
return op.distance();
241+
return detail::geos_safe_distance(geos_linestring_.get(), o.geos_polygon_.get());
244242
}
245243
template <typename T> template <typename U>
246244
double LineString<T>::distance(const Point<U>& o) const {
247-
geos::operation::distance::DistanceOp op(geos_linestring_.get(), o.geos_point_.get());
248-
return op.distance();
245+
return detail::geos_safe_distance(geos_linestring_.get(), o.geos_point_.get());
249246
}
250247
template <typename T> template <typename U>
251248
double LineString<T>::distance(const MultiLineString<U>& o) const {
252-
geos::operation::distance::DistanceOp op(geos_linestring_.get(), o.geos_mls_.get()); return op.distance();
249+
return detail::geos_safe_distance(geos_linestring_.get(), o.geos_mls_.get());
253250
}
254251

255252
// Python: shapely/geometry/base.py predicates L753-L813

‎shapely/geometry/point.h‎

Lines changed: 3 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -203,20 +203,17 @@ std::tuple<std::vector<T>, std::vector<T>> Point<T>::xy() const {
203203

204204
template <typename T>
205205
double Point<T>::distance(const Point& other) const {
206-
geos::operation::distance::DistanceOp op(geos_point_.get(), other.geos_point_.get());
207-
return op.distance();
206+
return detail::geos_safe_distance(geos_point_.get(), other.geos_point_.get());
208207
}
209208

210209
template <typename T> template <typename U>
211210
double Point<T>::distance(const LineString<U>& other) const {
212-
geos::operation::distance::DistanceOp op(geos_point_.get(), other.geos_linestring_.get());
213-
return op.distance();
211+
return detail::geos_safe_distance(geos_point_.get(), other.geos_linestring_.get());
214212
}
215213

216214
template <typename T> template <typename U>
217215
double Point<T>::distance(const Polygon<U>& other) const {
218-
geos::operation::distance::DistanceOp op(geos_point_.get(), other.geos_polygon_.get());
219-
return op.distance();
216+
return detail::geos_safe_distance(geos_point_.get(), other.geos_polygon_.get());
220217
}
221218

222219
// Python: shapely/geometry/base.py::buffer:L541

‎shapely/geometry/polygon.h‎

Lines changed: 8 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -458,18 +458,15 @@ void LinearRing<T>::normalize() { geos_ring_->normalize(); }
458458

459459
template <typename T>
460460
double LinearRing<T>::distance(const LinearRing& o) const {
461-
geos::operation::distance::DistanceOp op(geos_ring_.get(), o.geos_ring_.get());
462-
return op.distance();
461+
return detail::geos_safe_distance(geos_ring_.get(), o.geos_ring_.get());
463462
}
464463
template <typename T> template <typename U>
465464
double LinearRing<T>::distance(const Point<U>& o) const {
466-
geos::operation::distance::DistanceOp op(geos_ring_.get(), o.geos_point_.get());
467-
return op.distance();
465+
return detail::geos_safe_distance(geos_ring_.get(), o.geos_point_.get());
468466
}
469467
template <typename T> template <typename U>
470468
double LinearRing<T>::distance(const LineString<U>& o) const {
471-
geos::operation::distance::DistanceOp op(geos_ring_.get(), o.geos_linestring_.get());
472-
return op.distance();
469+
return detail::geos_safe_distance(geos_ring_.get(), o.geos_linestring_.get());
473470
}
474471

475472
// -- buffer ------------------------------------------------------------------
@@ -721,26 +718,23 @@ double Polygon<T>::area() const { return geos_polygon_->getArea(); }
721718

722719
template <typename T>
723720
double Polygon<T>::distance(const Polygon& o) const {
724-
geos::operation::distance::DistanceOp op(geos_polygon_.get(), o.geos_polygon_.get());
725-
return op.distance();
721+
return detail::geos_safe_distance(geos_polygon_.get(), o.geos_polygon_.get());
726722
}
727723
template <typename T> template <typename U>
728724
double Polygon<T>::distance(const LineString<U>& o) const {
729-
geos::operation::distance::DistanceOp op(geos_polygon_.get(), o.geos_linestring_.get());
730-
return op.distance();
725+
return detail::geos_safe_distance(geos_polygon_.get(), o.geos_linestring_.get());
731726
}
732727
template <typename T> template <typename U>
733728
double Polygon<T>::distance(const Point<U>& o) const {
734-
geos::operation::distance::DistanceOp op(geos_polygon_.get(), o.geos_point_.get());
735-
return op.distance();
729+
return detail::geos_safe_distance(geos_polygon_.get(), o.geos_point_.get());
736730
}
737731
template <typename T> template <typename U>
738732
double Polygon<T>::distance(const MultiLineString<U>& o) const {
739-
geos::operation::distance::DistanceOp op(geos_polygon_.get(), o.geos_mls_.get()); return op.distance();
733+
return detail::geos_safe_distance(geos_polygon_.get(), o.geos_mls_.get());
740734
}
741735
template <typename T> template <typename U>
742736
double Polygon<T>::distance(const MultiPolygon<U>& o) const {
743-
geos::operation::distance::DistanceOp op(geos_polygon_.get(), o.geos_mp_.get()); return op.distance();
737+
return detail::geos_safe_distance(geos_polygon_.get(), o.geos_mp_.get());
744738
}
745739

746740
// Python: shapely/geometry/base.py predicates L753-L813

‎tests/module.cpp‎

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -388,8 +388,9 @@ PYBIND11_MODULE(shapelycpp, m) {
388388
// Project / interpolate (from pycpp)
389389
// ======================================================================
390390
m.def("project_linestring_point", &project_ls_pt);
391-
m.def("interpolate_linestring", &interpolate_ls,
392-
py::arg("ls"), py::arg("distance"), py::arg("normalized") = false);
391+
m.def("interpolate_linestring", [](const LineString<double>& ls, double dist, bool normalized) {
392+
return interpolate_ls(ls, dist, normalized);
393+
}, py::arg("ls"), py::arg("distance"), py::arg("normalized") = false);
393394

394395
// ======================================================================
395396
// Intersection area (from pycpp)

‎tests/report.csv‎

Lines changed: 4 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -249,15 +249,7 @@ ops,nearest_points,unknown,f64,poly_ls_touch,PASS,0
249249
ops,nearest_points,unknown,f64,poly_ls_near,PASS,0
250250
ops,nearest_points_ls_pt,unknown,f64,ls_pt_above,PASS,0
251251
ops,nearest_points_ls_pt,unknown,f64,ls_pt_endpoint,PASS,0
252-
errors,point,unknown,f64,nan,PASS,0
253-
errors,point,unknown,f64,inf,PASS,0
254-
errors,linestring,unknown,f64,nan,PASS,0
255-
errors,linestring,unknown,f64,inf,PASS,0
256-
errors,linestring,unknown,f64,ninf,PASS,0
257-
errors,polygon,unknown,f64,nan,PASS,0
258-
errors,polygon,unknown,f64,inf,PASS,0
259-
errors,linearring,unknown,f64,nan,PASS,0
260-
errors,linearring,unknown,f64,inf,PASS,0
261-
errors,multipoint,unknown,f64,nan,PASS,0
262-
errors,multilinestring,unknown,f64,nan,PASS,0
263-
errors,multipolygon,unknown,f64,nan,PASS,0
252+
nan,linestring,unknown,f64,nan_ls,PASS,0
253+
nan,class.ls.length,unknown,f64,nan_len,PASS,0
254+
nan,class.ls.is_simple,unknown,f64,nan_simple,PASS,0
255+
nan,class.ls.is_valid,unknown,f64,nan_valid,PASS,0

0 commit comments

Comments
 (0)