fix: array_to_double_vec() hardcoded *2 corrupts (N,≥3) arrays; add C-contiguous guard

- array_to_double_vec: use row stride (r*cols) instead of flat rows*2 copy,
  correctly extracting only x,y columns from (N,≥3) arrays.
- Add ensure_c_contiguous() helper that checks buffer strides match
  C-order (row-major) layout — rejects non-contiguous views like
  arr[:,:2] with a clear ValueError directing caller to use .copy().
- Apply contiguity check to all 18 factory entry points (typed and
  generic overloads of point/linestring/polygon/linearring/multipoint/
  multilinestring/multipolygon).
- Fix point(py::array&) to handle both 1D [x,y] and 2D arrays directly
  instead of going through array_to_double_vec.
- tests: add TestArrayShape3D with 24 tests covering (N,3) extraction,
  row/col layout consistency, non-contiguous rejection, and .copy()
  workaround.

Author: peng.li24

pycpp/geometry_py.h

@@ -43,21 +43,57 @@ py::array_t<T> native_to_array(const T* data, size_t rows, size_t cols) {
     return result;
 }
 
+/// Throw if a py::buffer_info is not C-contiguous (row-major, no gaps).
+/// Protects against silent data corruption when callers pass sliced views
+/// (e.g. arr[:, :2], arr[::2], arr.T) whose memory strides don't match shape.
+inline void ensure_c_contiguous(const py::buffer_info& buf) {
+    if (buf.ndim == 0) return;
+    // Last dimension stride must equal element size
+    if (buf.strides[buf.ndim - 1] != buf.itemsize) {
+        throw std::invalid_argument(
+            "array must be C-contiguous; use .copy() on sliced/transposed views "
+            "(e.g. arr.copy() or arr[:, :2].copy())");
+    }
+    // Each higher dimension stride must equal next-dim-size × next-stride
+    for (ssize_t i = buf.ndim - 2; i >= 0; --i) {
+        ssize_t expected = buf.shape[i + 1] * buf.strides[i + 1];
+        if (buf.strides[i] != expected) {
+            throw std::invalid_argument(
+                "array must be C-contiguous; use .copy() on sliced/transposed views "
+                "(e.g. arr.copy() or arr[:, :2].copy())");
+        }
+    }
+}
+
 /// Convert any py::array coordinate buffer to std::vector<double>.
+/// Extracts the first 2 columns per row — safe for (N,2) and (N,≥3) arrays.
 /// dtype double → direct copy; float32/other → cast & widen to double.
 inline std::vector<double> array_to_double_vec(const py::array& arr) {
     auto buf = arr.request();
-    size_t sz = static_cast<size_t>(buf.shape[0]) * 2;
-    std::vector<double> tmp(sz);
+    if (buf.ndim != 2) {
+        throw std::invalid_argument("array_to_double_vec: expected 2D array");
+    }
+    ensure_c_contiguous(buf);
+    size_t rows = static_cast<size_t>(buf.shape[0]);
+    size_t cols = static_cast<size_t>(buf.shape[1]);
+    std::vector<double> tmp(rows * 2);
     if (arr.dtype().is(py::dtype::of<double>())) {
         const double* src = static_cast<const double*>(buf.ptr);
-        std::copy(src, src + sz, tmp.begin());
+        for (size_t r = 0; r < rows; ++r) {
+            size_t base = r * cols;
+            tmp[r * 2]     = src[base];
+            tmp[r * 2 + 1] = src[base + 1];
+        }
     } else {
         // float32 or other → cast through pybind11 to float32, then widen
         auto f32 = py::cast<py::array_t<float>>(arr);
         auto fbuf = f32.request();
         const float* src = static_cast<const float*>(fbuf.ptr);
-        for (size_t i = 0; i < sz; ++i) tmp[i] = static_cast<double>(src[i]);
+        for (size_t r = 0; r < rows; ++r) {
+            size_t base = r * cols;
+            tmp[r * 2]     = static_cast<double>(src[base]);
+            tmp[r * 2 + 1] = static_cast<double>(src[base + 1]);
+        }
     }
     return tmp;
 }
@@ -72,27 +108,43 @@ inline Point<float>  point(float x, float y)  { return Point<float>(x, y); }
 
 inline Point<double> point(const py::array_t<double>& arr) {
     auto buf = arr.request();
+    ensure_c_contiguous(buf);
     const double* p = static_cast<const double*>(buf.ptr);
     return Point<double>(buf.size > 0 ? p[0] : 0, buf.size > 1 ? p[1] : 0);
 }
 inline Point<float> point(const py::array_t<float>& arr) {
     auto buf = arr.request();
+    ensure_c_contiguous(buf);
     const float* p = static_cast<const float*>(buf.ptr);
     return Point<float>(buf.size > 0 ? p[0] : 0, buf.size > 1 ? p[1] : 0);
 }
 // auto-double: accept any dtype (int, unknown), always produce Point<double>
+// Handles 1D [x,y] and 2D (N,≥2) arrays — takes first 2 elements.
 inline Point<double> point(const py::array& arr) {
-    auto tmp = array_to_double_vec(arr);
-    return Point<double>(tmp[0], tmp[1]);
+    auto buf = arr.request();
+    ensure_c_contiguous(buf);
+    if (buf.size < 2) return Point<double>(0, 0);
+    if (arr.dtype().is(py::dtype::of<double>())) {
+        const double* p = static_cast<const double*>(buf.ptr);
+        return Point<double>(p[0], p[1]);
+    } else {
+        // float32 or int/other → cast through pybind11 to float32, then widen
+        auto f32 = py::cast<py::array_t<float>>(arr);
+        auto fbuf = f32.request();
+        const float* p = static_cast<const float*>(fbuf.ptr);
+        return Point<double>(static_cast<double>(p[0]), static_cast<double>(p[1]));
+    }
 }
 
 // -- LineString --
 inline LineString<double> linestring(const py::array_t<double>& arr) {
     auto buf = arr.request();
+    ensure_c_contiguous(buf);
     return LineString<double>(static_cast<const double*>(buf.ptr), buf.shape[0], buf.shape[1]);
 }
 inline LineString<float> linestring(const py::array_t<float>& arr) {
     auto buf = arr.request();
+    ensure_c_contiguous(buf);
     return LineString<float>(static_cast<const float*>(buf.ptr), buf.shape[0], buf.shape[1]);
 }
 inline LineString<double> linestring(const py::array& arr) {
@@ -104,10 +156,12 @@ inline LineString<double> linestring(const py::array& arr) {
 // -- Polygon --
 inline Polygon<double> polygon(const py::array_t<double>& arr) {
     auto buf = arr.request();
+    ensure_c_contiguous(buf);
     return Polygon<double>(static_cast<const double*>(buf.ptr), buf.shape[0], buf.shape[1]);
 }
 inline Polygon<float> polygon(const py::array_t<float>& arr) {
     auto buf = arr.request();
+    ensure_c_contiguous(buf);
     return Polygon<float>(static_cast<const float*>(buf.ptr), buf.shape[0], buf.shape[1]);
 }
 inline Polygon<double> polygon(const py::array& arr) {
@@ -119,6 +173,7 @@ inline Polygon<double> polygon(const py::array& arr) {
 // -- LinearRing (double only for now) --
 inline LinearRing<double> linearring(const py::array_t<double>& arr) {
     auto buf = arr.request();
+    ensure_c_contiguous(buf);
     return LinearRing<double>(static_cast<const double*>(buf.ptr), buf.shape[0], buf.shape[1]);
 }
 inline LinearRing<double> linearring(const py::array& arr) {
@@ -130,10 +185,12 @@ inline LinearRing<double> linearring(const py::array& arr) {
 // -- MultiPoint: single array of shape (n_pts, 2) --
 inline MultiPoint<double> multipoint(const py::array_t<double>& arr) {
     auto buf = arr.request();
+    ensure_c_contiguous(buf);
     return MultiPoint<double>(static_cast<const double*>(buf.ptr), buf.shape[0], buf.shape[1]);
 }
 inline MultiPoint<float> multipoint(const py::array_t<float>& arr) {
     auto buf = arr.request();
+    ensure_c_contiguous(buf);
     return MultiPoint<float>(static_cast<const float*>(buf.ptr), buf.shape[0], buf.shape[1]);
 }
 inline MultiPoint<double> multipoint(const py::array& arr) {
@@ -147,6 +204,7 @@ inline MultiLineString<double> multilinestring(const std::vector<py::array_t<dou
     MultiLineString<double> mls;
     for (auto& arr : arrays) {
         auto buf = arr.request();
+        ensure_c_contiguous(buf);
         mls.add_line(static_cast<const double*>(buf.ptr), buf.shape[0], buf.shape[1]);
     }
     return mls;
@@ -155,6 +213,7 @@ inline MultiLineString<float> multilinestring(const std::vector<py::array_t<floa
     MultiLineString<float> mls;
     for (auto& arr : arrays) {
         auto buf = arr.request();
+        ensure_c_contiguous(buf);
         mls.add_line(static_cast<const float*>(buf.ptr), buf.shape[0], buf.shape[1]);
     }
     return mls;
@@ -174,6 +233,7 @@ inline MultiPolygon<double> multipolygon(const std::vector<py::array_t<double>>&
     MultiPolygon<double> mp;
     for (auto& arr : arrays) {
         auto buf = arr.request();
+        ensure_c_contiguous(buf);
         mp.add_polygon(static_cast<const double*>(buf.ptr), buf.shape[0], buf.shape[1]);
     }
     return mp;
@@ -182,6 +242,7 @@ inline MultiPolygon<float> multipolygon(const std::vector<py::array_t<float>>& a
     MultiPolygon<float> mp;
     for (auto& arr : arrays) {
         auto buf = arr.request();
+        ensure_c_contiguous(buf);
         mp.add_polygon(static_cast<const float*>(buf.ptr), buf.shape[0], buf.shape[1]);
     }
     return mp;

tests/test_api.py

@@ -1050,3 +1050,180 @@ def test_random_ls_poly_f64(self, cpp):
             i_cpp = cpp.intersects_linestring_polygon(cpp.linestring(l), cpp.polygon(c))
             i_py = PyLineString(l).intersects(PyPolygon(c))
             assert i_cpp == i_py
+
+
+# ==============================================================================
+# 3-COLUMN ARRAY SHAPE HANDLING — verify (N,≥3) → x,y extraction
+# ==============================================================================
+# Regression tests for the array_to_double_vec() fix (hardcoded rows*2 → stride).
+# Covers: generic py::array overload (int dtype), typed double overload, point,
+# polygon, multipoint; plus row/column layout consistency and non-contiguous view guard.
+
+class TestArrayShape3D:
+    """Verify (N,≥3) arrays correctly extract x,y columns across all factory overloads."""
+
+    # -- (N,2) regression: generic overload (int32) still works ----------------------
+
+    def test_2col_generic(self, cpp):
+        coords = np.array([[0, 1], [2, 3], [4, 5]], dtype=np.int32)
+        ls = cpp.linestring(coords)
+        assert ls.coords() == [(0.0, 1.0), (2.0, 3.0), (4.0, 5.0)]
+
+    # -- (N,3) generic overload: z column correctly ignored -------------------------
+
+    def test_3col_generic(self, cpp):
+        coords = np.array([[10, 20, 999], [30, 40, 888], [50, 60, 777]], dtype=np.int32)
+        ls = cpp.linestring(coords)
+        assert ls.coords() == [(10.0, 20.0), (30.0, 40.0), (50.0, 60.0)]
+
+    # -- (N,3) typed double overload: stride by shape[1]=3 --------------------------
+
+    def test_3col_typed_double(self, cpp):
+        coords = np.array([[10.0, 20.0, 99.0], [30.0, 40.0, 88.0], [50.0, 60.0, 77.0]],
+                          dtype=np.float64)
+        ls = cpp.linestring(coords)
+        assert ls.coords() == [(10.0, 20.0), (30.0, 40.0), (50.0, 60.0)]
+
+    # -- Equality: (N,3) produces identical coords to (N,2) -------------------------
+
+    def test_3col_equals_2col_generic(self, cpp):
+        coords_3 = np.array([[1, 2, 99], [3, 4, 88], [5, 6, 77]], dtype=np.int32)
+        coords_2 = np.array([[1.0, 2.0], [3.0, 4.0], [5.0, 6.0]], dtype=np.float64)
+        ls3 = cpp.linestring(coords_3)
+        ls2 = cpp.linestring(coords_2)
+        assert ls3.coords() == ls2.coords()
+        # Distance to self should also match
+        assert ls3.length() == ls2.length()
+
+    def test_3col_equals_2col_typed(self, cpp):
+        coords_3 = np.array([[1.0, 2.0, 99.0], [3.0, 4.0, 88.0], [5.0, 6.0, 77.0]],
+                            dtype=np.float64)
+        coords_2 = np.array([[1.0, 2.0], [3.0, 4.0], [5.0, 6.0]], dtype=np.float64)
+        ls3 = cpp.linestring(coords_3)
+        ls2 = cpp.linestring(coords_2)
+        assert ls3.coords() == ls2.coords()
+        assert ls3.length() == ls2.length()
+
+    # -- Polygon (N,3) generic ------------------------------------------------------
+
+    def test_polygon_3col_generic(self, cpp):
+        coords = np.array([[0, 0, 1], [10, 0, 2], [10, 10, 3], [0, 10, 4], [0, 0, 5]],
+                          dtype=np.int32)
+        poly = cpp.polygon(coords)
+        ext = cpp.polygon_exterior(poly)
+        # First 5 points must match (GEOS may auto-close to 6 points)
+        assert ext.tolist()[:5] == [[0.0, 0.0], [10.0, 0.0], [10.0, 10.0], [0.0, 10.0], [0.0, 0.0]]
+
+    def test_polygon_3col_typed(self, cpp):
+        coords = np.array([[0.0, 0.0, 1.0], [10.0, 0.0, 2.0], [10.0, 10.0, 3.0],
+                           [0.0, 10.0, 4.0], [0.0, 0.0, 5.0]], dtype=np.float64)
+        poly = cpp.polygon(coords)
+        ext = cpp.polygon_exterior(poly)
+        assert ext.tolist()[:5] == [[0.0, 0.0], [10.0, 0.0], [10.0, 10.0], [0.0, 10.0], [0.0, 0.0]]
+
+    # -- Point: 1D [x,y] and 2D (N,≥2) int arrays (auto-double overload) -------------
+
+    def test_point_1d_int(self, cpp):
+        p = cpp.point(np.array([3, 4], dtype=np.int32))
+        assert p.coords() == [(3.0, 4.0)]
+
+    def test_point_2d_int(self, cpp):
+        p = cpp.point(np.array([[5, 6]], dtype=np.int32))
+        assert p.coords() == [(5.0, 6.0)]
+
+    def test_point_2d_3col_int(self, cpp):
+        """Point from (N,3) int: takes first 2 elements (row 0, col 0-1)."""
+        p = cpp.point(np.array([[7, 8, 999], [9, 10, 888]], dtype=np.int32))
+        assert p.coords() == [(7.0, 8.0)]
+
+    # -- MultiPoint (N,3) generic ---------------------------------------------------
+
+    def test_multipoint_3col_generic(self, cpp):
+        mp = cpp.multipoint(np.array([[1, 2, 99], [4, 6, 88]], dtype=np.int32))
+        # Verify via distance: (1,2) should be a point in the multipoint
+        p = cpp.point(1.0, 2.0)
+        assert mp.distance(p) == 0.0
+
+    # -- LinearRing (N,3) -----------------------------------------------------------
+
+    def test_linearring_3col(self, cpp, C):
+        ring = C.linearring(np.array([[0, 0, 1], [1, 0, 2], [1, 1, 3], [0, 0, 4]],
+                                      dtype=np.float64))
+        assert ring.is_closed() == True
+        assert ring.is_ring() == True
+
+    # -- Row/column layout consistency ----------------------------------------------
+
+    def test_row_col_layout(self, cpp):
+        """Verify Python arr[i,j] ↔ C++ ptr[i*shape[1]+j] mapping."""
+        # (3,4) array: row-major [1,2,3,4, 5,6,7,8, 9,10,11,12]
+        arr = np.array([[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12]], dtype=np.float64)
+        ls = cpp.linestring(arr)
+        # C++ takes coords_[i*4+0], coords_[i*4+1] → first 2 columns per row
+        assert ls.coords() == [(1.0, 2.0), (5.0, 6.0), (9.0, 10.0)]
+
+        # (3,3) direct: stride=3, skip z
+        arr3 = np.array([[10, 20, 99], [30, 40, 88], [50, 60, 77]], dtype=np.float64)
+        ls3 = cpp.linestring(arr3)
+        assert ls3.coords() == [(10.0, 20.0), (30.0, 40.0), (50.0, 60.0)]
+
+    # -- poly_exterior output shape --------------------------------------------------
+
+    def test_exterior_shape(self, cpp):
+        coords = np.array([[0, 0], [10, 0], [10, 10], [0, 10]], dtype=np.float64)
+        poly = cpp.polygon(coords)
+        ext = cpp.polygon_exterior(poly)
+        assert ext.ndim == 2
+        assert ext.shape[1] == 2  # always 2 columns
+        assert ext.dtype == np.float64
+
+    # -- Non-contiguous view rejection (fail-fast with clear error) ----------------
+
+    def test_noncontiguous_view_rejected_typed(self, cpp):
+        """Typed overloads throw on non-contiguous views (e.g. [:, :2])."""
+        arr = np.array([[10, 20, 99], [30, 40, 88], [50, 60, 77]], dtype=np.float64)
+        view = arr[:, :2]
+        assert not view.flags['C_CONTIGUOUS']
+        with pytest.raises(ValueError, match="C-contiguous"):
+            cpp.linestring(view)
+
+    def test_noncontiguous_view_rejected_generic(self, cpp):
+        """Generic (int dtype) overloads throw on non-contiguous views."""
+        arr = np.array([[10, 20, 99], [30, 40, 88], [50, 60, 77]], dtype=np.int32)
+        view = arr[:, :2]
+        assert not view.flags['C_CONTIGUOUS']
+        with pytest.raises(ValueError, match="C-contiguous"):
+            cpp.linestring(view)
+
+    def test_noncontiguous_copy_works_typed(self, cpp):
+        """.copy() on a non-contiguous view produces correct results."""
+        arr = np.array([[10, 20, 99], [30, 40, 88], [50, 60, 77]], dtype=np.float64)
+        cpy = arr[:, :2].copy()
+        assert cpy.flags['C_CONTIGUOUS']
+        ls = cpp.linestring(cpy)
+        assert ls.coords() == [(10.0, 20.0), (30.0, 40.0), (50.0, 60.0)]
+
+    def test_noncontiguous_copy_works_generic(self, cpp):
+        """.copy() on a non-contiguous view works for generic overload."""
+        arr = np.array([[10, 20, 99], [30, 40, 88], [50, 60, 77]], dtype=np.int32)
+        cpy = arr[:, :2].copy()
+        assert cpy.flags['C_CONTIGUOUS']
+        ls = cpp.linestring(cpy)
+        assert ls.coords() == [(10.0, 20.0), (30.0, 40.0), (50.0, 60.0)]
+
+    @pytest.mark.parametrize("factory, arr_dtype", [
+        ("point", np.float64),
+        ("point", np.int32),
+        ("polygon", np.float64),
+        ("polygon", np.int32),
+        ("multipoint", np.float64),
+        ("multipoint", np.int32),
+    ])
+    def test_noncontiguous_rejected_all_factories(self, cpp, factory, arr_dtype):
+        """All factory functions reject non-contiguous views."""
+        arr = np.array([[10, 20, 99], [30, 40, 88], [50, 60, 77]], dtype=arr_dtype)
+        view = arr[:, :2]
+        assert not view.flags['C_CONTIGUOUS']
+        fn = getattr(cpp, factory)
+        with pytest.raises(ValueError, match="C-contiguous"):
+            fn(view)