ianthomas23 · Aug 28, 2022 · Aug 31, 2022 · Aug 31, 2022 · Sep 1, 2022 · Sep 2, 2022
diff --git a/lib/matplotlib/tests/test_triangulation.py b/lib/matplotlib/tests/test_triangulation.py
@@ -1167,43 +1167,43 @@ def test_internal_cpp_api():
     # C++ Triangulation.
     with pytest.raises(
             TypeError,
-            match=r'function takes exactly 7 arguments \(0 given\)'):
+            match=r'__init__\(\): incompatible constructor arguments.'):
         mpl._tri.Triangulation()
 
     with pytest.raises(
             ValueError, match=r'x and y must be 1D arrays of the same length'):
-        mpl._tri.Triangulation([], [1], [[]], None, None, None, False)
+        mpl._tri.Triangulation([], [1], [[]], (), (), (), False)
 
     x = [0, 1, 1]
     y = [0, 0, 1]
     with pytest.raises(
             ValueError,
             match=r'triangles must be a 2D array of shape \(\?,3\)'):
-        mpl._tri.Triangulation(x, y, [[0, 1]], None, None, None, False)
+        mpl._tri.Triangulation(x, y, [[0, 1]], (), (), (), False)
 
     tris = [[0, 1, 2]]
     with pytest.raises(
             ValueError,
             match=r'mask must be a 1D array with the same length as the '
                   r'triangles array'):
-        mpl._tri.Triangulation(x, y, tris, [0, 1], None, None, False)
+        mpl._tri.Triangulation(x, y, tris, [0, 1], (), (), False)
 
     with pytest.raises(
             ValueError, match=r'edges must be a 2D array with shape \(\?,2\)'):
-        mpl._tri.Triangulation(x, y, tris, None, [[1]], None, False)
+        mpl._tri.Triangulation(x, y, tris, (), [[1]], (), False)
 
     with pytest.raises(
             ValueError,
             match=r'neighbors must be a 2D array with the same shape as the '
                   r'triangles array'):
-        mpl._tri.Triangulation(x, y, tris, None, None, [[-1]], False)
+        mpl._tri.Triangulation(x, y, tris, (), (), [[-1]], False)
 
-    triang = mpl._tri.Triangulation(x, y, tris, None, None, None, False)
+    triang = mpl._tri.Triangulation(x, y, tris, (), (), (), False)
 
     with pytest.raises(
             ValueError,
-            match=r'z array must have same length as triangulation x and y '
-                  r'array'):
+            match=r'z must be a 1D array with the same length as the '
+                  r'triangulation x and y arrays'):
         triang.calculate_plane_coefficients([])
 
     with pytest.raises(
@@ -1215,7 +1215,7 @@ def test_internal_cpp_api():
     # C++ TriContourGenerator.
     with pytest.raises(
             TypeError,
-            match=r'function takes exactly 2 arguments \(0 given\)'):
+            match=r'__init__\(\): incompatible constructor arguments.'):
         mpl._tri.TriContourGenerator()
 
     with pytest.raises(
@@ -1233,7 +1233,8 @@ def test_internal_cpp_api():
 
     # C++ TrapezoidMapTriFinder.
     with pytest.raises(
-            TypeError, match=r'function takes exactly 1 argument \(0 given\)'):
+            TypeError,
+            match=r'__init__\(\): incompatible constructor arguments.'):
         mpl._tri.TrapezoidMapTriFinder()
 
     trifinder = mpl._tri.TrapezoidMapTriFinder(triang)

diff --git a/lib/matplotlib/tri/triangulation.py b/lib/matplotlib/tri/triangulation.py
@@ -120,8 +120,12 @@ def get_cpp_triangulation(self):
         from matplotlib import _tri
         if self._cpp_triangulation is None:
             self._cpp_triangulation = _tri.Triangulation(
-                self.x, self.y, self.triangles, self.mask, self._edges,
-                self._neighbors, not self.is_delaunay)
+                # For unset arrays use empty tuple which has size of zero.
+                self.x, self.y, self.triangles,
+                self.mask if self.mask is not None else (),
+                self._edges if self._edges is not None else (),
+                self._neighbors if self._neighbors is not None else (),
+                not self.is_delaunay)
         return self._cpp_triangulation
 
     def get_masked_triangles(self):

diff --git a/setup.py b/setup.py
@@ -302,6 +302,8 @@ def make_release_tree(self, base_dir, files):
     setup_requires=[
         "certifi>=2020.06.20",
         "numpy>=1.19",
+        "pybind11>=2.6",
+        "setuptools>=42",
         "setuptools_scm>=7",
     ],
     install_requires=[

diff --git a/setupext.py b/setupext.py
@@ -16,6 +16,7 @@
 import textwrap
 import urllib.request
 
+from pybind11.setup_helpers import Pybind11Extension
 from setuptools import Distribution, Extension
 
 _log = logging.getLogger(__name__)
@@ -432,10 +433,10 @@ def get_extensions(self):
         add_libagg_flags_and_sources(ext)
         yield ext
         # qhull
-        ext = Extension(
+        ext = Pybind11Extension(
             "matplotlib._qhull", ["src/_qhull_wrapper.cpp"],
+            cxx_std=11,
             define_macros=[("MPL_DEVNULL", os.devnull)])
-        add_numpy_flags(ext)
         Qhull.add_flags(ext)
         yield ext
         # tkagg
@@ -452,12 +453,12 @@ def get_extensions(self):
         add_libagg_flags(ext)
         yield ext
         # tri
-        ext = Extension(
+        ext = Pybind11Extension(
             "matplotlib._tri", [
                 "src/tri/_tri.cpp",
                 "src/tri/_tri_wrapper.cpp",
-            ])
-        add_numpy_flags(ext)
+            ],
+            cxx_std=11)
         yield ext
         # ttconv
         ext = Extension(

diff --git a/src/_qhull_wrapper.cpp b/src/_qhull_wrapper.cpp
@@ -5,9 +5,9 @@
  * triangulation, construct an instance of the matplotlib.tri.Triangulation
  * class without specifying a triangles array.
  */
-#define PY_SSIZE_T_CLEAN
-#include "Python.h"
-#include "numpy_cpp.h"
+#include <pybind11/pybind11.h>
+#include <pybind11/numpy.h>
+
 #ifdef _MSC_VER
 /* The Qhull header does not declare this as extern "C", but only MSVC seems to
  * do name mangling on global variables. We thus need to declare this before
@@ -16,18 +16,28 @@ extern "C" {
 extern const char qh_version[];
 }
 #endif
+
 #include "libqhull_r/qhull_ra.h"
 #include <cstdio>
+#include <sstream>
 #include <vector>
 
-
 #ifndef MPL_DEVNULL
 #error "MPL_DEVNULL must be defined as the OS-equivalent of /dev/null"
 #endif
 
 #define STRINGIFY(x) STR(x)
 #define STR(x) #x
 
+namespace py = pybind11;
+
+// Input numpy array class.
+typedef py::array_t<double, py::array::c_style | py::array::forcecast> CoordArray;
+
+// Output numpy array class.
+typedef py::array_t<int> IndexArray;
+
+
 
 static const char* qhull_error_msg[6] = {
     "",                     /* 0 = qh_ERRnone */
@@ -64,11 +74,11 @@ get_facet_neighbours(const facetT* facet, std::vector<int>& tri_indices,
 /* Return true if the specified points arrays contain at least 3 unique points,
  * or false otherwise. */
 static bool
-at_least_3_unique_points(npy_intp npoints, const double* x, const double* y)
+at_least_3_unique_points(py::ssize_t npoints, const double* x, const double* y)
 {
-    int i;
-    const int unique1 = 0;  /* First unique point has index 0. */
-    int unique2 = 0;        /* Second unique point index is 0 until set. */
+    py::ssize_t i;
+    const py::ssize_t unique1 = 0;  /* First unique point has index 0. */
+    py::ssize_t unique2 = 0;        /* Second unique point index is 0 until set. */
 
     if (npoints < 3) {
         return false;
@@ -125,8 +135,8 @@ class QhullInfo {
 /* Delaunay implementation method.
  * If hide_qhull_errors is true then qhull error messages are discarded;
  * if it is false then they are written to stderr. */
-static PyObject*
-delaunay_impl(npy_intp npoints, const double* x, const double* y,
+py::tuple
+delaunay_impl(py::ssize_t npoints, const double* x, const double* y,
               bool hide_qhull_errors)
 {
     qhT qh_qh;                  /* qh variable type and name must be like */
@@ -179,11 +189,12 @@ delaunay_impl(npy_intp npoints, const double* x, const double* y,
     exitcode = qh_new_qhull(qh, ndim, (int)npoints, points.data(), False,
                             (char*)"qhull d Qt Qbb Qc Qz", NULL, error_file);
     if (exitcode != qh_ERRnone) {
-        PyErr_Format(PyExc_RuntimeError,
-                     "Error in qhull Delaunay triangulation calculation: %s (exitcode=%d)%s",
-                     qhull_error_msg[exitcode], exitcode,
-                     hide_qhull_errors ? "; use python verbose option (-v) to see original qhull error." : "");
-        return NULL;
+        std::ostringstream oss;
+        oss << "Error in qhull Delaunay triangulation calculation: " <<
+            qhull_error_msg[exitcode] << " (exitcode=" << exitcode << ")";
+        if (hide_qhull_errors)
+            oss << "; use python verbose option (-v) to see original qhull error.";
+        throw std::runtime_error(oss.str());
     }
 
     /* Split facets so that they only have 3 points each. */
@@ -204,12 +215,12 @@ delaunay_impl(npy_intp npoints, const double* x, const double* y,
     std::vector<int> tri_indices(max_facet_id+1);
 
     /* Allocate Python arrays to return. */
-    npy_intp dims[2] = {ntri, 3};
-    numpy::array_view<int, ndim> triangles(dims);
-    int* triangles_ptr = triangles.data();
+    int dims[2] = {ntri, 3};
+    IndexArray triangles(dims);
+    int* triangles_ptr = triangles.mutable_data();
 
-    numpy::array_view<int, ndim> neighbors(dims);
-    int* neighbors_ptr = neighbors.data();
+    IndexArray neighbors(dims);
+    int* neighbors_ptr = neighbors.mutable_data();
 
     /* Determine triangles array and set tri_indices array. */
     i = 0;
@@ -238,103 +249,56 @@ delaunay_impl(npy_intp npoints, const double* x, const double* y,
         }
     }
 
-    PyObject* tuple = PyTuple_New(2);
-    if (tuple == 0) {
-        throw std::runtime_error("Failed to create Python tuple");
-    }
-
-    PyTuple_SET_ITEM(tuple, 0, triangles.pyobj());
-    PyTuple_SET_ITEM(tuple, 1, neighbors.pyobj());
-    return tuple;
+    return py::make_tuple(triangles, neighbors);
 }
 
 /* Process Python arguments and call Delaunay implementation method. */
-static PyObject*
-delaunay(PyObject *self, PyObject *args)
+py::tuple
+delaunay(const CoordArray& x, const CoordArray& y)
 {
-    numpy::array_view<double, 1> xarray;
-    numpy::array_view<double, 1> yarray;
-    PyObject* ret;
-    npy_intp npoints;
-    const double* x;
-    const double* y;
-
-    if (!PyArg_ParseTuple(args, "O&O&",
-                          &xarray.converter_contiguous, &xarray,
-                          &yarray.converter_contiguous, &yarray)) {
-        return NULL;
-    }
+    if (x.ndim() != 1 || y.ndim() != 1)
+        throw std::invalid_argument("x and y must be 1D arrays");
 
-    npoints = xarray.dim(0);
-    if (npoints != yarray.dim(0)) {
-        PyErr_SetString(PyExc_ValueError,
-                        "x and y must be 1D arrays of the same length");
-        return NULL;
-    }
+    auto npoints = x.shape(0);
+    if (y.shape(0) != npoints)
+        throw std::invalid_argument("x and y must be 1D arrays of the same length");
 
-    if (npoints < 3) {
-        PyErr_SetString(PyExc_ValueError,
-                        "x and y arrays must have a length of at least 3");
-        return NULL;
-    }
-
-    x = xarray.data();
-    y = yarray.data();
-
-    if (!at_least_3_unique_points(npoints, x, y)) {
-        PyErr_SetString(PyExc_ValueError,
-                        "x and y arrays must consist of at least 3 unique points");
-        return NULL;
-    }
+    if (npoints < 3)
+        throw std::invalid_argument("x and y arrays must have a length of at least 3");
 
-    CALL_CPP("qhull.delaunay",
-             (ret = delaunay_impl(npoints, x, y, Py_VerboseFlag == 0)));
+    if (!at_least_3_unique_points(npoints, x.data(), y.data()))
+        throw std::invalid_argument("x and y arrays must consist of at least 3 unique points");
 
-    return ret;
+    return delaunay_impl(npoints, x.data(), y.data(), Py_VerboseFlag == 0);
 }
 
 /* Return qhull version string for assistance in debugging. */
-static PyObject*
-version(PyObject *self, PyObject *arg)
+std::string
+version()
 {
-    return PyBytes_FromString(qh_version);
+    return std::string(qh_version);
 }
 
-static PyMethodDef qhull_methods[] = {
-    {"delaunay", delaunay, METH_VARARGS,
-     "delaunay(x, y, /)\n"
-     "--\n\n"
-     "Compute a Delaunay triangulation.\n"
-     "\n"
-     "Parameters\n"
-     "----------\n"
-     "x, y : 1d arrays\n"
-     "    The coordinates of the point set, which must consist of at least\n"
-     "    three unique points.\n"
-     "\n"
-     "Returns\n"
-     "-------\n"
-     "triangles, neighbors : int arrays, shape (ntri, 3)\n"
-     "    Indices of triangle vertices and indices of triangle neighbors.\n"
-    },
-    {"version", version, METH_NOARGS,
-     "version()\n--\n\n"
-     "Return the qhull version string."},
-    {NULL, NULL, 0, NULL}
-};
-
-static struct PyModuleDef qhull_module = {
-    PyModuleDef_HEAD_INIT,
-    "qhull", "Computing Delaunay triangulations.\n", -1, qhull_methods
-};
-
-#pragma GCC visibility push(default)
-
-PyMODINIT_FUNC
-PyInit__qhull(void)
-{
-    import_array();
-    return PyModule_Create(&qhull_module);
+PYBIND11_MODULE(_qhull, m) {
+    m.doc() = "Computing Delaunay triangulations.\n";
+
+    m.def("delaunay", &delaunay,
+        "delaunay(x, y, /)\n"
+        "--\n\n"
+        "Compute a Delaunay triangulation.\n"
+        "\n"
+        "Parameters\n"
+        "----------\n"
+        "x, y : 1d arrays\n"
+        "    The coordinates of the point set, which must consist of at least\n"
+        "    three unique points.\n"
+        "\n"
+        "Returns\n"
+        "-------\n"
+        "triangles, neighbors : int arrays, shape (ntri, 3)\n"
+        "    Indices of triangle vertices and indices of triangle neighbors.\n");
+
+    m.def("version", &version,
+        "version()\n--\n\n"
+        "Return the qhull version string.");
 }
-
-#pragma GCC visibility pop