5
#include "VoronoiDiagramGenerator.h"
6
#include "delaunay_utils.h"
7
#include "natneighbors.h"
8
#include "numpy/noprefix.h"
14
static void reorder_edges(int npoints, int ntriangles,
16
int *edge_db, int *tri_edges, int *tri_nbrs)
18
int neighbors[3], nodes[3];
22
for (i=0; i<ntriangles; i++) {
23
nodes[0] = INDEX2(edge_db, INDEX3(tri_edges,i,0), 0);
24
nodes[1] = INDEX2(edge_db, INDEX3(tri_edges,i,0), 1);
25
tmp = INDEX2(edge_db, INDEX3(tri_edges,i,1), 0);
26
if (tmp == nodes[0]) {
28
nodes[2] = INDEX2(edge_db, INDEX3(tri_edges,i,1), 1);
29
} else if (tmp == nodes[1]) {
31
nodes[2] = INDEX2(edge_db, INDEX3(tri_edges,i,1), 1);
32
} else if (INDEX2(edge_db, INDEX3(tri_edges,i,1), 1) == nodes[0]) {
40
if (ONRIGHT(x[nodes[0]], y[nodes[0]],
41
x[nodes[1]], y[nodes[1]],
42
x[nodes[2]], y[nodes[2]])) {
43
// flip to make counter-clockwise
50
// I worked it out on paper. You're just gonna have to trust me on this.
51
if (!case1 && !case2) {
52
neighbors[0] = INDEX3(tri_nbrs, i, 1);
53
neighbors[1] = INDEX3(tri_nbrs, i, 2);
54
neighbors[2] = INDEX3(tri_nbrs, i, 0);
55
} else if (case1 && !case2) {
56
neighbors[0] = INDEX3(tri_nbrs, i, 2);
57
neighbors[1] = INDEX3(tri_nbrs, i, 1);
58
neighbors[2] = INDEX3(tri_nbrs, i, 0);
59
} else if (!case1 && case2) {
60
neighbors[0] = INDEX3(tri_nbrs, i, 1);
61
neighbors[1] = INDEX3(tri_nbrs, i, 0);
62
neighbors[2] = INDEX3(tri_nbrs, i, 2);
64
neighbors[0] = INDEX3(tri_nbrs, i, 2);
65
neighbors[1] = INDEX3(tri_nbrs, i, 0);
66
neighbors[2] = INDEX3(tri_nbrs, i, 1);
69
// Not trusting me? Okay, let's go through it:
70
// We have three edges to deal with and three nodes. Without loss
71
// of generality, let's label the nodes A, B, and C with (A, B)
72
// forming the first edge in the order they arrive on input.
73
// Then there are eight possibilities as to how the other edge-tuples
74
// may be labeled, but only two variations that are going to affect the
81
// The distinction is whether A is in the second edge or B is.
82
// This is the test "case1" above.
84
// The second test we need to perform is for counter-clockwiseness.
85
// Again, there are only two variations that will affect the outcome:
86
// either ABC is counter-clockwise, or it isn't. In the former case,
87
// we're done setting the node order, we just need to associate the
88
// appropriate neighbor triangles with their opposite nodes, something
89
// which can be done by inspection. In the latter case, to order the
90
// nodes counter-clockwise, we only have to switch B and C to get
91
// nodes ACB. Then we simply set the neighbor list by inspection again.
101
// ABC ACB -+- node order
104
INDEX3(tri_edges,i,0) = nodes[0];
105
INDEX3(tri_edges,i,1) = nodes[1];
106
INDEX3(tri_edges,i,2) = nodes[2];
107
INDEX3(tri_nbrs,i,0) = neighbors[0];
108
INDEX3(tri_nbrs,i,1) = neighbors[1];
109
INDEX3(tri_nbrs,i,2) = neighbors[2];
113
static PyObject* getMesh(int npoints, double *x, double *y)
115
PyObject *vertices, *edge_db, *tri_edges, *tri_nbrs;
117
int tri0, tri1, reg0, reg1;
118
double tri0x, tri0y, tri1x, tri1y;
119
int length, numtri, i, j;
121
int *edge_db_ptr, *tri_edges_ptr, *tri_nbrs_ptr;
122
double *vertices_ptr;
123
VoronoiDiagramGenerator vdg;
125
vdg.generateVoronoi(x, y, npoints, -100, 100, -100, 100, 0);
126
vdg.getNumbers(length, numtri);
130
edge_db = PyArray_SimpleNew(2, dim, PyArray_INT);
131
if (!edge_db) goto fail;
132
edge_db_ptr = (int*)PyArray_DATA(edge_db);
135
vertices = PyArray_SimpleNew(2, dim, PyArray_DOUBLE);
136
if (!vertices) goto fail;
137
vertices_ptr = (double*)PyArray_DATA(vertices);
140
tri_edges = PyArray_SimpleNew(2, dim, PyArray_INT);
141
if (!tri_edges) goto fail;
142
tri_edges_ptr = (int*)PyArray_DATA(tri_edges);
144
tri_nbrs = PyArray_SimpleNew(2, dim, PyArray_INT);
145
if (!tri_nbrs) goto fail;
146
tri_nbrs_ptr = (int*)PyArray_DATA(tri_nbrs);
148
for (i=0; i<(3*numtri); i++) {
149
tri_edges_ptr[i] = tri_nbrs_ptr[i] = -1;
152
vdg.resetEdgeListIter();
154
while (vdg.getNextDelaunay(tri0, tri0x, tri0y, tri1, tri1x, tri1y, reg0, reg1)) {
156
INDEX2(edge_db_ptr,i,0) = reg0;
157
INDEX2(edge_db_ptr,i,1) = reg1;
159
INDEX2(vertices_ptr,tri0,0) = tri0x;
160
INDEX2(vertices_ptr,tri0,1) = tri0y;
161
for (j=0; j<3; j++) {
162
if (INDEX3(tri_edges_ptr,tri0,j) == i) break;
163
if (INDEX3(tri_edges_ptr,tri0,j) == -1) {
164
INDEX3(tri_edges_ptr,tri0,j) = i;
165
INDEX3(tri_nbrs_ptr,tri0,j) = tri1;
171
INDEX2(vertices_ptr,tri1,0) = tri1x;
172
INDEX2(vertices_ptr,tri1,1) = tri1y;
173
for (j=0; j<3; j++) {
174
if (INDEX3(tri_edges_ptr,tri1,j) == i) break;
175
if (INDEX3(tri_edges_ptr,tri1,j) == -1) {
176
INDEX3(tri_edges_ptr,tri1,j) = i;
177
INDEX3(tri_nbrs_ptr,tri1,j) = tri0;
184
// tri_edges contains lists of edges; convert to lists of nodes in
185
// counterclockwise order and reorder tri_nbrs to match. Each node
186
// corresponds to the edge opposite it in the triangle.
187
reorder_edges(npoints, numtri, x, y, edge_db_ptr, tri_edges_ptr,
190
temp = Py_BuildValue("(OOOO)", vertices, edge_db, tri_edges, tri_nbrs);
191
if (!temp) goto fail;
195
Py_DECREF(tri_edges);
201
Py_XDECREF(vertices);
203
Py_XDECREF(tri_edges);
204
Py_XDECREF(tri_nbrs);
208
static PyObject *linear_planes(int ntriangles, double *x, double *y, double *z,
215
double x02, y02, z02, x12, y12, z12, xy0212;
217
dims[0] = ntriangles;
219
planes = PyArray_SimpleNew(2, dims, PyArray_DOUBLE);
220
if (!planes) return NULL;
221
planes_ptr = (double *)PyArray_DATA(planes);
223
for (i=0; i<ntriangles; i++) {
224
x02 = x[INDEX3(nodes,i,0)] - x[INDEX3(nodes,i,2)];
225
y02 = y[INDEX3(nodes,i,0)] - y[INDEX3(nodes,i,2)];
226
z02 = z[INDEX3(nodes,i,0)] - z[INDEX3(nodes,i,2)];
227
x12 = x[INDEX3(nodes,i,1)] - x[INDEX3(nodes,i,2)];
228
y12 = y[INDEX3(nodes,i,1)] - y[INDEX3(nodes,i,2)];
229
z12 = z[INDEX3(nodes,i,1)] - z[INDEX3(nodes,i,2)];
233
INDEX3(planes_ptr,i,0) = (z02 - z12 * xy0212) / (x02 - x12 * xy0212);
234
INDEX3(planes_ptr,i,1) = (z12 - INDEX3(planes_ptr,i,0)*x12) / y12;
235
INDEX3(planes_ptr,i,2) = (z[INDEX3(nodes,i,2)] -
236
INDEX3(planes_ptr,i,0)*x[INDEX3(nodes,i,2)] -
237
INDEX3(planes_ptr,i,1)*y[INDEX3(nodes,i,2)]);
240
INDEX3(planes_ptr,i,1) = (z02 - z12 * xy0212) / (y02 - y12 * xy0212);
241
INDEX3(planes_ptr,i,0) = (z12 - INDEX3(planes_ptr,i,1)*y12) / x12;
242
INDEX3(planes_ptr,i,2) = (z[INDEX3(nodes,i,2)] -
243
INDEX3(planes_ptr,i,0)*x[INDEX3(nodes,i,2)] -
244
INDEX3(planes_ptr,i,1)*y[INDEX3(nodes,i,2)]);
248
return (PyObject*)planes;
251
static double linear_interpolate_single(double targetx, double targety,
252
double *x, double *y, int *nodes, int *neighbors,
253
PyObject *planes, double defvalue, int start_triangle, int *end_triangle)
256
planes_ptr = (double*)PyArray_DATA(planes);
257
if (start_triangle == -1) start_triangle = 0;
258
*end_triangle = walking_triangles(start_triangle, targetx, targety,
259
x, y, nodes, neighbors);
260
if (*end_triangle == -1) return defvalue;
261
return (targetx*INDEX3(planes_ptr,*end_triangle,0) +
262
targety*INDEX3(planes_ptr,*end_triangle,1) +
263
INDEX3(planes_ptr,*end_triangle,2));
266
static PyObject *linear_interpolate_grid(double x0, double x1, int xsteps,
267
double y0, double y1, int ysteps,
268
PyObject *planes, double defvalue,
269
int npoints, double *x, double *y, int *nodes, int *neighbors)
272
double dx, dy, targetx, targety;
273
int rowtri, coltri, tri;
280
z = PyArray_SimpleNew(2, dims, PyArray_DOUBLE);
282
z_ptr = (double*)PyArray_DATA(z);
284
dx = (x1 - x0) / (xsteps-1);
285
dy = (y1 - y0) / (ysteps-1);
288
for (iy=0; iy<ysteps; iy++) {
289
targety = y0 + dy*iy;
290
rowtri = walking_triangles(rowtri, x0, targety, x, y, nodes, neighbors);
292
for (ix=0; ix<xsteps; ix++) {
293
targetx = x0 + dx*ix;
294
INDEXN(z_ptr, xsteps, iy, ix) = linear_interpolate_single(
296
x, y, nodes, neighbors, planes, defvalue, tri, &coltri);
297
if (coltri != -1) tri = coltri;
304
static PyObject *compute_planes_method(PyObject *self, PyObject *args)
306
PyObject *pyx, *pyy, *pyz, *pynodes;
307
PyArrayObject *x, *y, *z, *nodes;
308
int npoints, ntriangles;
312
if (!PyArg_ParseTuple(args, "OOOO", &pyx, &pyy, &pyz, &pynodes)) {
315
x = (PyArrayObject*)PyArray_ContiguousFromObject(pyx, PyArray_DOUBLE, 1, 1);
316
y = (PyArrayObject*)PyArray_ContiguousFromObject(pyy, PyArray_DOUBLE, 1, 1);
317
z = (PyArrayObject*)PyArray_ContiguousFromObject(pyz, PyArray_DOUBLE, 1, 1);
318
if ((!x || !y) || !z) {
319
PyErr_SetString(PyExc_ValueError, "x,y,z must be 1-D arrays of floats");
322
npoints = PyArray_DIM(x, 0);
323
if ((PyArray_DIM(y, 0) != npoints) || (PyArray_DIM(z, 0) != npoints)) {
324
PyErr_SetString(PyExc_ValueError, "x,y,z arrays must be of equal length");
327
nodes = (PyArrayObject*)PyArray_ContiguousFromObject(pynodes, PyArray_INT, 2, 2);
329
PyErr_SetString(PyExc_ValueError, "nodes must be a 2-D array of ints");
332
ntriangles = PyArray_DIM(nodes, 0);
333
if (PyArray_DIM(nodes, 1) != 3) {
334
PyErr_SetString(PyExc_ValueError, "nodes must have shape (ntriangles, 3)");
338
planes = linear_planes(ntriangles, (double*)PyArray_DATA(x),
339
(double*)PyArray_DATA(y), (double*)PyArray_DATA(z), (int*)PyArray_DATA(nodes));
356
static PyObject *linear_interpolate_method(PyObject *self, PyObject *args)
358
double x0, x1, y0, y1, defvalue;
360
PyObject *pyplanes, *pyx, *pyy, *pynodes, *pyneighbors, *grid;
363
PyArrayObject *planes, *x, *y, *nodes, *neighbors;
365
if (!PyArg_ParseTuple(args, "ddiddidOOOOO", &x0, &x1, &xsteps, &y0, &y1, &ysteps,
366
&defvalue, &pyplanes, &pyx, &pyy, &pynodes, &pyneighbors)) {
369
x = (PyArrayObject*)PyArray_ContiguousFromObject(pyx, PyArray_DOUBLE, 1, 1);
370
y = (PyArrayObject*)PyArray_ContiguousFromObject(pyy, PyArray_DOUBLE, 1, 1);
372
PyErr_SetString(PyExc_ValueError, "x,y must be 1-D arrays of floats");
375
npoints = PyArray_DIM(x, 0);
376
if (PyArray_DIM(y, 0) != npoints) {
377
PyErr_SetString(PyExc_ValueError, "x,y arrays must be of equal length");
380
planes = (PyArrayObject*)PyArray_ContiguousFromObject(pyplanes, PyArray_DOUBLE, 2, 2);
382
PyErr_SetString(PyExc_ValueError, "planes must be a 2-D array of floats");
385
nodes = (PyArrayObject*)PyArray_ContiguousFromObject(pynodes, PyArray_INT, 2, 2);
387
PyErr_SetString(PyExc_ValueError, "nodes must be a 2-D array of ints");
390
neighbors = (PyArrayObject*)PyArray_ContiguousFromObject(pyneighbors, PyArray_INT, 2, 2);
392
PyErr_SetString(PyExc_ValueError, "neighbors must be a 2-D array of ints");
396
grid = linear_interpolate_grid(x0, x1, xsteps, y0, y1, ysteps,
397
(PyObject*)planes, defvalue, npoints,
398
(double*)PyArray_DATA(x), (double*)PyArray_DATA(y),
399
(int*)PyArray_DATA(nodes), (int*)PyArray_DATA(neighbors));
405
Py_DECREF(neighbors);
414
Py_XDECREF(neighbors);
424
Py_XDECREF(centers);\
426
Py_XDECREF(neighbors);\
430
static PyObject *nn_interpolate_unstructured_method(PyObject *self, PyObject *args)
432
PyObject *pyx, *pyy, *pyz, *pycenters, *pynodes, *pyneighbors, *pyintx, *pyinty;
433
PyArrayObject *x, *y, *z, *centers, *nodes, *neighbors, *intx, *inty, *intz;
435
int size, npoints, ntriangles;
437
if (!PyArg_ParseTuple(args, "OOdOOOOOO", &pyintx, &pyinty, &defvalue,
438
&pyx, &pyy, &pyz, &pycenters, &pynodes, &pyneighbors)) {
441
x = (PyArrayObject*)PyArray_ContiguousFromObject(pyx, PyArray_DOUBLE, 1, 1);
442
y = (PyArrayObject*)PyArray_ContiguousFromObject(pyy, PyArray_DOUBLE, 1, 1);
443
z = (PyArrayObject*)PyArray_ContiguousFromObject(pyz, PyArray_DOUBLE, 1, 1);
444
if ((!x || !y) || !z) {
445
PyErr_SetString(PyExc_ValueError, "x,y must be 1-D arrays of floats");
449
npoints = PyArray_DIM(x, 0);
450
if ((PyArray_DIM(y, 0) != npoints) || (PyArray_DIM(z, 0) != npoints)) {
451
PyErr_SetString(PyExc_ValueError, "x,y,z arrays must be of equal length");
455
centers = (PyArrayObject*)PyArray_ContiguousFromObject(pycenters, PyArray_DOUBLE, 2, 2);
457
PyErr_SetString(PyExc_ValueError, "centers must be a 2-D array of ints");
461
nodes = (PyArrayObject*)PyArray_ContiguousFromObject(pynodes, PyArray_INT, 2, 2);
463
PyErr_SetString(PyExc_ValueError, "nodes must be a 2-D array of ints");
467
neighbors = (PyArrayObject*)PyArray_ContiguousFromObject(pyneighbors, PyArray_INT, 2, 2);
469
PyErr_SetString(PyExc_ValueError, "neighbors must be a 2-D array of ints");
473
ntriangles = PyArray_DIM(neighbors, 0);
474
if ((PyArray_DIM(nodes, 0) != ntriangles) ||
475
(PyArray_DIM(centers, 0) != ntriangles)) {
476
PyErr_SetString(PyExc_ValueError, "centers,nodes,neighbors must be of equal length");
480
intx = (PyArrayObject*)PyArray_ContiguousFromObject(pyintx, PyArray_DOUBLE, 0, 0);
481
inty = (PyArrayObject*)PyArray_ContiguousFromObject(pyinty, PyArray_DOUBLE, 0, 0);
482
if (!intx || !inty) {
483
PyErr_SetString(PyExc_ValueError, "intx,inty must be arrays of floats");
486
if (intx->nd != inty->nd) {
487
PyErr_SetString(PyExc_ValueError, "intx,inty must have same shapes");
490
for (int i=0; i<intx->nd; i++) {
491
if (intx->dimensions[i] != inty->dimensions[i]) {
492
PyErr_SetString(PyExc_ValueError, "intx,inty must have same shapes");
496
intz = (PyArrayObject*)PyArray_SimpleNew(intx->nd, intx->dimensions, PyArray_DOUBLE);
499
NaturalNeighbors nn(npoints, ntriangles,
500
(double*)PyArray_DATA(x), (double*)PyArray_DATA(y),
501
(double*)PyArray_DATA(centers), (int*)PyArray_DATA(nodes),
502
(int*)PyArray_DATA(neighbors));
503
size = PyArray_Size((PyObject*)intx);
504
nn.interpolate_unstructured((double*)PyArray_DATA(z), size,
505
(double*)PyArray_DATA(intx), (double*)PyArray_DATA(inty),
506
(double*)PyArray_DATA(intz), defvalue);
508
return (PyObject*)intz;
517
Py_XDECREF(centers);\
519
Py_XDECREF(neighbors);\
522
static PyObject *nn_interpolate_method(PyObject *self, PyObject *args)
524
PyObject *pyx, *pyy, *pyz, *pycenters, *pynodes, *pyneighbors, *grid;
525
PyArrayObject *x, *y, *z, *centers, *nodes, *neighbors;
526
double x0, x1, y0, y1, defvalue;
528
int npoints, ntriangles;
531
if (!PyArg_ParseTuple(args, "ddiddidOOOOOO", &x0, &x1, &xsteps,
532
&y0, &y1, &ysteps, &defvalue, &pyx, &pyy, &pyz, &pycenters, &pynodes,
536
x = (PyArrayObject*)PyArray_ContiguousFromObject(pyx, PyArray_DOUBLE, 1, 1);
537
y = (PyArrayObject*)PyArray_ContiguousFromObject(pyy, PyArray_DOUBLE, 1, 1);
538
z = (PyArrayObject*)PyArray_ContiguousFromObject(pyz, PyArray_DOUBLE, 1, 1);
539
if ((!x || !y) || !z) {
540
PyErr_SetString(PyExc_ValueError, "x,y must be 1-D arrays of floats");
544
npoints = PyArray_DIM(x, 0);
545
if (PyArray_DIM(y, 0) != npoints) {
546
PyErr_SetString(PyExc_ValueError, "x,y arrays must be of equal length");
550
centers = (PyArrayObject*)PyArray_ContiguousFromObject(pycenters, PyArray_DOUBLE, 2, 2);
552
PyErr_SetString(PyExc_ValueError, "centers must be a 2-D array of ints");
556
nodes = (PyArrayObject*)PyArray_ContiguousFromObject(pynodes, PyArray_INT, 2, 2);
558
PyErr_SetString(PyExc_ValueError, "nodes must be a 2-D array of ints");
562
neighbors = (PyArrayObject*)PyArray_ContiguousFromObject(pyneighbors, PyArray_INT, 2, 2);
564
PyErr_SetString(PyExc_ValueError, "neighbors must be a 2-D array of ints");
568
ntriangles = PyArray_DIM(neighbors, 0);
569
if ((PyArray_DIM(nodes, 0) != ntriangles) ||
570
(PyArray_DIM(centers, 0) != ntriangles)) {
571
PyErr_SetString(PyExc_ValueError, "centers,nodes,neighbors must be of equal length");
575
//goto succeed; // XXX: Can't cross NaturalNeighbors instantiation with goto
579
grid = PyArray_SimpleNew(2, dims, PyArray_DOUBLE);
580
if (!grid) {CLEANUP} // goto fail;
582
NaturalNeighbors nn(npoints, ntriangles,
583
(double*)PyArray_DATA(x), (double*)PyArray_DATA(y),
584
(double*)PyArray_DATA(centers), (int*)PyArray_DATA(nodes),
585
(int*)PyArray_DATA(neighbors));
586
nn.interpolate_grid((double*)PyArray_DATA(z),
589
(double*)PyArray_DATA(grid),
597
Py_DECREF(neighbors);
604
static PyObject *delaunay_method(PyObject *self, PyObject *args)
606
PyObject *pyx, *pyy, *mesh;
607
PyArrayObject *x, *y;
610
if (!PyArg_ParseTuple(args, "OO", &pyx, &pyy)) {
613
x = (PyArrayObject*)PyArray_ContiguousFromObject(pyx, PyArray_DOUBLE, 1, 1);
614
y = (PyArrayObject*)PyArray_ContiguousFromObject(pyy, PyArray_DOUBLE, 1, 1);
616
PyErr_SetString(PyExc_ValueError, "x,y must be 1-D arrays of floats");
620
npoints = PyArray_DIM(x, 0);
621
if (PyArray_DIM(y, 0) != npoints) {
622
PyErr_SetString(PyExc_ValueError, "x and y must have the same length");
626
mesh = getMesh(npoints, (double*)(x->data), (double*)(y->data));
628
if (!mesh) goto fail;
641
static PyMethodDef delaunay_methods[] = {
642
{"delaunay", (PyCFunction)delaunay_method, METH_VARARGS,
643
"Compute the Delaunay triangulation of a cloud of 2-D points.\n\n"
644
"circumcenters, edges, tri_points, tri_neighbors = delaunay(x, y)\n\n"
645
"x, y -- shape-(npoints,) arrays of floats giving the X and Y coordinates of the points\n"
646
"circumcenters -- shape-(numtri,2) array of floats giving the coordinates of the\n"
647
" circumcenters of each triangle (numtri being the number of triangles)\n"
648
"edges -- shape-(nedges,2) array of integers giving the indices into x and y\n"
649
" of each edge in the triangulation\n"
650
"tri_points -- shape-(numtri,3) array of integers giving the indices into x and y\n"
651
" of each node in each triangle\n"
652
"tri_neighbors -- shape-(numtri,3) array of integers giving the indices into circumcenters\n"
653
" tri_points, and tri_neighbors of the neighbors of each triangle\n"},
654
{"compute_planes", (PyCFunction)compute_planes_method, METH_VARARGS,
656
{"linear_interpolate_grid", (PyCFunction)linear_interpolate_method, METH_VARARGS,
658
{"nn_interpolate_grid", (PyCFunction)nn_interpolate_method, METH_VARARGS,
660
{"nn_interpolate_unstructured", (PyCFunction)nn_interpolate_unstructured_method, METH_VARARGS,
662
{NULL, NULL, 0, NULL}
666
PyMODINIT_FUNC init_delaunay(void)
669
m = Py_InitModule3("_delaunay", delaunay_methods,
670
"Tools for computing the Delaunay triangulation and some operations on it.\n"