BUG: Delaunay/qhull segfault on self-pass

[tylerjereddy]

I noticed this segfault on main while investigating gh-21279. python dev.py test -- -k "test_find_simplex_self_segfault". It is clearly an abuse of the call signature but we generally should intercept before segfaulting.

--- a/scipy/spatial/tests/test_qhull.py
+++ b/scipy/spatial/tests/test_qhull.py
@@ -194,6 +194,12 @@ class TestUtilities:
             j = qhull.tsearch(tri, p[:2])
             assert_equal(i, j)
 
+    def test_find_simplex_self_segfault(self):
+        # segfaults before bug fix
+        points2d = np.array([[0, 0], [0, 1.1], [1, 0], [1, 1]])
+        tri = qhull.Delaunay(points2d)
+        tri.find_simplex(tri)
+
     def test_plane_distance(self):
         # Compare plane distance from hyperplane equations obtained from Qhull
         # to manually computed plane equations

[ev-br]

To clarify, the abuse of signature you refer to is that find_simplex expects an array and segfaults if given non-array arguments?

Looking at https://github.com/scipy/scipy/blob/main/scipy/spatial/_qhull.pyx#L2113, it does something strange indeed: it calls np.asanyarray, but are subclasses supported really. So I guess in your call, xi becomes an object array and things go askew from this point on?

Also tri.find_simplex(1) also segfaults for me locally, even though I'd expect the dimensionality check to ValueError instead.

[tylerjereddy]

There seem to be some things to address yeah. To be honest, I was kind of following the docs since our current devdocs docstring for find_simplex has the triangulation object/info as a valid parameter (see below screen capture as well). So I thought "that's weird" and tried it.

[ev-br]

Huh, indeed. But find_simplex does not take the tri parameter, so the docs are actively misleading :-).