#define PY_SSIZE_T_CLEAN

#include "Python.h"

#include "numpy_cpp.h"

#ifdef _MSC_VER

extern "C" {

extern const char qh_version[];

}

#endif

#include "libqhull_r/qhull_ra.h"

#include <cstdio>

#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

static const char* qhull_error_msg[6] = {

"", /* 0 = qh_ERRnone */

"input inconsistency", /* 1 = qh_ERRinput */

"singular input data", /* 2 = qh_ERRsingular */

"precision error", /* 3 = qh_ERRprec */

"insufficient memory", /* 4 = qh_ERRmem */

"internal error"}; /* 5 = qh_ERRqhull */

static void

get_facet_vertices(qhT* qh, const facetT* facet, int indices[3])

}

static void

get_facet_neighbours(const facetT* facet, std::vector<int>& tri_indices,

}

static bool

at_least_3_unique_points(npy_intp npoints, const double* x, const double* y)

}

class QhullInfo {

};

static PyObject*

delaunay_impl(npy_intp npoints, const double* x, const double* y,

bool hide_qhull_errors)

{

qhT qh_qh; /* qh variable type and name must be like */

qhT* qh = &qh_qh; /* this for Qhull macros to work correctly. */

facetT* facet;

int i, ntri, max_facet_id;

int exitcode; /* Value returned from qh_new_qhull(). */

const int ndim = 2;

double x_mean = 0.0;

double y_mean = 0.0;

QHULL_LIB_CHECK

/* Allocate points. */

std::vector<coordT> points(npoints * ndim);

/* Determine mean x, y coordinates. */

for (i = 0; i < npoints; ++i) {

x_mean += x[i];

y_mean += y[i];

}

x_mean /= npoints;

y_mean /= npoints;

/* Prepare points array to pass to qhull. */

for (i = 0; i < npoints; ++i) {

points[2*i ] = x[i] - x_mean;

points[2*i+1] = y[i] - y_mean;

}

/* qhull expects a FILE* to write errors to. */

FILE* error_file = NULL;

if (hide_qhull_errors) {

/* qhull errors are ignored by writing to OS-equivalent of /dev/null.

error_file = fopen(STRINGIFY(MPL_DEVNULL), "w");

if (error_file == NULL) {

throw std::runtime_error("Could not open devnull");

}

}

else {

/* qhull errors written to stderr. */

error_file = stderr;

}

/* Perform Delaunay triangulation. */

QhullInfo info(error_file, qh);

qh_zero(qh, error_file);

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;

}

/* Split facets so that they only have 3 points each. */

qh_triangulate(qh);

/* Determine ntri and max_facet_id.

ntri = 0;

FORALLfacets {

if (!facet->upperdelaunay) {

++ntri;

}

}

max_facet_id = qh->facet_id - 1;

/* Create array to map facet id to triangle index. */

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();

numpy::array_view<int, ndim> neighbors(dims);

int* neighbors_ptr = neighbors.data();

/* Determine triangles array and set tri_indices array. */

i = 0;

FORALLfacets {

if (!facet->upperdelaunay) {

int indices[3];

tri_indices[facet->id] = i++;

get_facet_vertices(qh, facet, indices);

*triangles_ptr++ = (facet->toporient ? indices[0] : indices[2]);

*triangles_ptr++ = indices[1];

*triangles_ptr++ = (facet->toporient ? indices[2] : indices[0]);

}

else {

tri_indices[facet->id] = -1;

}

}

/* Determine neighbors array. */

FORALLfacets {

if (!facet->upperdelaunay) {

int indices[3];

get_facet_neighbours(facet, tri_indices, indices);

*neighbors_ptr++ = (facet->toporient ? indices[2] : indices[0]);

*neighbors_ptr++ = (facet->toporient ? indices[0] : indices[2]);

*neighbors_ptr++ = indices[1];

}

}

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;

}

static PyObject*

delaunay(PyObject *self, PyObject *args)

{

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;

}

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;

}

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;

}

CALL_CPP("qhull.delaunay",

(ret = delaunay_impl(npoints, x, y, Py_VerboseFlag == 0)));

return ret;

}

static PyObject*

version(PyObject *self, PyObject *arg)

{

return PyBytes_FromString(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)

}

#pragma GCC visibility pop