Create triangulation.c

philippeitis · web-flow · commit 2c4f18690e19 · 2019-12-15T01:28:33.000-08:00
triangulation.c provides fast implementations for some of the functions in learner/triangulation.py
diff --git a/adaptive/learner/triangulation.c b/adaptive/learner/triangulation.c
@@ -0,0 +1,272 @@
+#include "Python.h"
+#include "numpy/arrayobject.h"
+#include <math.h>
+#include <stdio.h>
+#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
+
+static PyObject *
+fast_norm(PyObject *self, PyObject *args)
+{
+    PyObject* vect;
+	if (!PyArg_ParseTuple(args, "O", &vect))
+		return NULL;
+
+    if (PyList_Check(vect)) {
+        Py_ssize_t numElements = PyObject_Length(vect);
+        if (numElements <= 0) {
+            PyErr_SetString(PyExc_ValueError, "fast_norm requires a list of length 1 or greater.");
+            return NULL;
+        }
+        double sum = 0;
+        for (Py_ssize_t i = 0; i < numElements; ++i) {
+            double val = PyFloat_AsDouble(PyList_GetItem(vect, i));
+            if (PyErr_Occurred()) {
+                return NULL;
+            }
+            sum += val * val;
+        }
+        return Py_BuildValue("f", sqrt(sum));
+    }
+
+    if (PyTuple_Check(vect)) {
+        Py_ssize_t numElements = PyTuple_Size(vect);
+        if (numElements <= 0) {
+            PyErr_SetString(PyExc_ValueError, "fast_norm requires a list of length 1 or greater.");
+            return NULL;
+        }
+        double sum = 0;
+        for (Py_ssize_t i = 0; i < numElements; ++i) {
+            double val = PyFloat_AsDouble(PyTuple_GetItem(vect, i));
+            if (PyErr_Occurred()) {
+                return NULL;
+            }
+            sum += val * val;
+        }
+        return Py_BuildValue("f", sqrt(sum));
+    }
+
+    if (PyArray_Check(vect)) {
+        PyObject *npArray = PyArray_FROM_OTF(vect, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY);
+        if (npArray == NULL)
+            return NULL;
+        int nd = PyArray_NDIM(npArray);
+        if (nd != 1)
+            return NULL;
+        npy_intp * shape = PyArray_DIMS(npArray);
+        Py_ssize_t numElements = PyArray_SIZE(npArray);
+        double* data = (double *) PyArray_DATA(npArray);
+        double sum = 0;
+        for (Py_ssize_t i = 0; i < numElements; ++i) {
+            double val = data[i];
+            sum += val * val;
+        }
+        return Py_BuildValue("f", sqrt(sum));
+    }
+
+}
+
+static int get_tuple_elems(PyObject * tuple, double * first, double * second) {
+    PyObject * elem = PyTuple_GetItem(tuple, 0);
+    if (elem == NULL)
+        return 0;
+    *first = PyFloat_AsDouble(elem);
+    elem = PyTuple_GetItem(tuple, 1);
+    if (elem == NULL)
+        return 0;
+    *second = PyFloat_AsDouble(elem);
+    return 1;
+}
+
+static int get_triple_tuple_elems(PyObject * tuple, double * first, double * second, double* third) {
+    PyObject * elem = PyTuple_GetItem(tuple, 0);
+    if (elem == NULL)
+        return 0;
+    *first = PyFloat_AsDouble(elem);
+    elem = PyTuple_GetItem(tuple, 1);
+    if (elem == NULL)
+        return 0;
+    *second = PyFloat_AsDouble(elem);
+    elem = PyTuple_GetItem(tuple, 2);
+    if (elem == NULL)
+        return 0;
+    *third = PyFloat_AsDouble(elem);
+    return 1;
+}
+
+static PyObject *
+fast_2d_point_in_simplex(PyObject *self, PyObject *args, PyObject *kwargs)
+{
+    PyObject* point, * simplex;
+    double eps = 0.00000001;
+
+    // SystemError: bad format char passed to Py_BuildValue ??
+	static char *kwlist[] = {"point", "simplex", "eps", NULL};
+
+	if (!PyArg_ParseTupleAndKeywords(args, kwargs, "OO|d", kwlist, &point, &simplex, &eps))
+		return NULL;
+
+
+    if (!PyTuple_Check(point))
+        return NULL;
+
+    if (!PyList_Check(simplex))
+        return NULL;
+
+    double px, py, p0x, p0y, p1x, p1y, p2x, p2y;
+    if (!get_tuple_elems(point, &px, &py))
+        return NULL;
+
+    point = PyList_GetItem(simplex, 0);
+    if (point == NULL)
+        return NULL;
+
+    if (!get_tuple_elems(point, &p0x, &p0y))
+        return NULL;
+
+    point = PyList_GetItem(simplex, 1);
+    if (point == NULL)
+        return NULL;
+
+    if (!get_tuple_elems(point, &p1x, &p1y))
+        return NULL;
+
+    point = PyList_GetItem(simplex, 2);
+    if (point == NULL)
+        return NULL;
+
+    if (!get_tuple_elems(point, &p2x, &p2y))
+        return NULL;
+
+    double area = 0.5 * (-p1y * p2x + p0y * (p2x - p1x) + p1x * p2y + p0x * (p1y - p2y));
+    double s = 1 / (2 * area) * (p0y * p2x + (p2y - p0y) * px - p0x * p2y + (p0x - p2x) * py);
+    if ((s < -eps) || (s > 1 + eps))
+        return Py_BuildValue("O", Py_False);
+    double t = 1 / (2 * area) * (p0x * p1y + (p0y - p1y) * px - p0y * p1x + (p1x - p0x) * py);
+    return Py_BuildValue("O", (t >= -eps) && (s + t <= 1 + eps) ? Py_True : Py_False);
+}
+
+static PyObject *
+fast_2d_circumcircle(PyObject *self, PyObject *args)
+{
+    PyObject* points;
+	if (!PyArg_ParseTuple(args, "O", &points))
+		return NULL;
+
+    if (PyList_Check(points)) {
+        Py_ssize_t numElements = PyObject_Length(points);
+        if (numElements <= 0) {
+            PyErr_SetString(PyExc_ValueError, "fast_2d_circumcircle requires a list of length 1 or greater.");
+            return NULL;
+        }
+        double x0, y0, x1, y1, x2, y2;
+        if (!get_tuple_elems(PyList_GetItem(points, 0), &x0, &y0))
+            return NULL;
+        if (!get_tuple_elems(PyList_GetItem(points, 1), &x1, &y1))
+            return NULL;
+        if (!get_tuple_elems(PyList_GetItem(points, 2), &x2, &y2))
+            return NULL;
+
+        x1 -= x0;
+        y1 -= y0;
+
+        x2 -= x0;
+        y2 -= y0;
+
+        double l1 = x1 * x1 + y1 * y1;
+        double l2 = x2 * x2 + y2 * y2;
+        double dx = l1 * y2 - l2 * y1;
+        double dy = -l1 * x2 + l2 * x1;
+
+        double aa = 2 * (x1 * y2 - x2 * y1);
+        double x = dx / aa;
+        double y = dy / aa;
+        return Py_BuildValue("(ff)f", x + x0, y + y0, sqrt(x * x + y * y));
+    }
+    return NULL;
+}
+
+static PyObject *
+fast_3d_circumcircle(PyObject *self, PyObject *args)
+{
+    PyObject* points;
+	if (!PyArg_ParseTuple(args, "O", &points))
+		return NULL;
+
+    if (PyList_Check(points)) {
+        Py_ssize_t numElements = PyObject_Length(points);
+        if (numElements <= 0) {
+            PyErr_SetString(PyExc_ValueError, "fast_2d_circumcircle requires a list of length 1 or greater.");
+            return NULL;
+        }
+        double x0, y0, z0, x1, y1, z1, x2, y2, z2, x3, y3, z3;
+        if (!get_triple_tuple_elems(PyList_GetItem(points, 0), &x0, &y0, &z0))
+            return NULL;
+        if (!get_triple_tuple_elems(PyList_GetItem(points, 1), &x1, &y1, &z1))
+            return NULL;
+        if (!get_triple_tuple_elems(PyList_GetItem(points, 2), &x2, &y2, &z2))
+            return NULL;
+        if (!get_triple_tuple_elems(PyList_GetItem(points, 3), &x3, &y3, &z3))
+            return NULL;
+
+        x1 -= x0;
+        y1 -= y0;
+        z1 -= z0;
+
+        x2 -= x0;
+        y2 -= y0;
+        z2 -= z0;
+
+        x3 -= x0;
+        y3 -= y0;
+        z3 -= z0;
+
+        double l1 = x1 * x1 + y1 * y1 + z1 * z1;
+        double l2 = x2 * x2 + y2 * y2 + z2 * z2;
+        double l3 = x3 * x3 + y3 * y3 + z3 * z3;
+        double dx = l1 * (y2 * z3 - z2 * y3) - l2 * (y1 * z3 - z1 * y3) + l3 * (y1 * z2 - z1 * y2);
+        double dy = -l1 * (x2 * z3 - z2 * x3) - l2 * (x1 * z3 - z1 * x3) + l3 * (x1 * z2 - z1 * x2);
+        double dz = l1 * (x2 * y3 - y2 * x3) - l2 * (x1 * y3 - y1 * x3) + l3 * (x1 * y2 - y1 * x2);
+
+        double aa = 2 * (x1 * (y2 * z3 - z2 * y3) - x2 * (y1 * z3 - z1 * y3) + x3 * (y1 * z2 - z1 * y2));
+        double x = dx / aa;
+        double y = dy / aa;
+        double z = dz / aa;
+        return Py_BuildValue("(fff)f", x + x0, y + y0, z + z0, sqrt(x * x + y * y + z * z));
+    }
+    return NULL;
+}
+
+
+static PyMethodDef triangulation_functions[] = {
+    {"fast_norm", fast_norm, METH_VARARGS,
+		"Returns the norm of the given array."},
+    {"fast_2d_circumcircle", fast_2d_circumcircle, METH_VARARGS,
+		"Returns the norm of the given array."},
+    {"fast_3d_circumcircle", fast_3d_circumcircle, METH_VARARGS,
+		"Returns the norm of the given array."},
+    {"fast_2d_point_in_simplex", (PyCFunction) fast_2d_point_in_simplex, METH_VARARGS | METH_KEYWORDS,
+		"Refer to docs."},
+	{NULL}
+};
+
+PyDoc_STRVAR(module_doc, "Fast implementations of triangulation functions.");
+
+static struct PyModuleDef moduledef = {
+	PyModuleDef_HEAD_INIT,
+	"triangulation",
+	module_doc,
+	-1,                 /* m_size */
+	triangulation_functions,   /* m_methods */
+	NULL,               /* m_reload (unused) */
+	NULL,               /* m_traverse */
+	NULL,               /* m_clear */
+	NULL                /* m_free */
+};
+
+PyObject *
+PyInit_triangulation(void)
+{
+    if(PyArray_API == NULL)
+        import_array();
+    return PyModule_Create(&moduledef);
+}
