← Back to branch summary

~ubuntu-branches/ubuntu/saucy/python-scipy/saucy

« back to all changes in this revision

Viewing changes to scipy/sandbox/delaunay/_delaunay.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Ondrej Certik
  • Date: 2008-06-16 22:58:01 UTC
  • mfrom: (2.1.24 intrepid)
  • Revision ID: james.westby@ubuntu.com-20080616225801-irdhrpcwiocfbcmt
Tags: 0.6.0-12
http://bugs.debian.org/489149
* The description updated to match the current SciPy (Closes: #489149).
* Standards-Version bumped to 3.8.0 (no action needed)
* Build-Depends: netcdf-dev changed to libnetcdf-dev

Show diffs side-by-side

added added

removed removed

Lines of Context:
scipy/sandbox/delaunay/_delaunay.cpp
 
1
#include "Python.h"
 
2
#include <stdlib.h>
 
3
#include <map>
 
4
 
 
5
#include "VoronoiDiagramGenerator.h"
 
6
#include "delaunay_utils.h"
 
7
#include "natneighbors.h"
 
8
#include "numpy/noprefix.h"
 
9
 
 
10
using namespace std;
 
11
 
 
12
extern "C" {
 
13
 
 
14
static void reorder_edges(int npoints, int ntriangles, 
 
15
    double *x, double *y, 
 
16
    int *edge_db, int *tri_edges, int *tri_nbrs)
 
17
{
 
18
    int neighbors[3], nodes[3];
 
19
    int i, tmp;
 
20
    int case1, case2;
 
21
 
 
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]) {
 
27
            case1 = 1;
 
28
            nodes[2] = INDEX2(edge_db, INDEX3(tri_edges,i,1), 1);
 
29
        } else if (tmp == nodes[1]) {
 
30
            case1 = 0;
 
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]) {
 
33
            case1 = 1;
 
34
            nodes[2] = tmp;
 
35
        } else {
 
36
            case1 = 0;
 
37
            nodes[2] = tmp;
 
38
        }
 
39
 
 
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
 
44
            tmp = nodes[2];
 
45
            nodes[2] = nodes[1];
 
46
            nodes[1] = tmp;
 
47
            case2 = 1;
 
48
        } else case2 = 0;
 
49
 
 
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);
 
63
        } else {
 
64
            neighbors[0] = INDEX3(tri_nbrs, i, 2);
 
65
            neighbors[1] = INDEX3(tri_nbrs, i, 0);
 
66
            neighbors[2] = INDEX3(tri_nbrs, i, 1);
 
67
        }
 
68
 
 
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
 
75
        // output:
 
76
        //
 
77
        // AB         AB
 
78
        // BC (CB)    AC (CA)
 
79
        // CA (AC)    BC (CB)
 
80
        //
 
81
        // The distinction is whether A is in the second edge or B is.
 
82
        // This is the test "case1" above.
 
83
        //
 
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.
 
92
        //
 
93
        //          CCW     CW
 
94
        // AB
 
95
        // BC       120     102 -+
 
96
        // CA                    |
 
97
        //                       +- neighbor order
 
98
        // AB                    |
 
99
        // AC       210     201 -+
 
100
        // BC
 
101
        //          ABC     ACB -+- node order
 
102
 
 
103
 
 
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];
 
110
    }
 
111
}
 
112
 
 
113
static PyObject* getMesh(int npoints, double *x, double *y)
 
114
{
 
115
    PyObject *vertices, *edge_db, *tri_edges, *tri_nbrs;
 
116
    PyObject *temp;
 
117
    int tri0, tri1, reg0, reg1;
 
118
    double tri0x, tri0y, tri1x, tri1y;
 
119
    int length, numtri, i, j;
 
120
    intp dim[MAX_DIMS];
 
121
    int *edge_db_ptr, *tri_edges_ptr, *tri_nbrs_ptr;
 
122
    double *vertices_ptr;
 
123
    VoronoiDiagramGenerator vdg;
 
124
 
 
125
    vdg.generateVoronoi(x, y, npoints, -100, 100, -100, 100, 0);
 
126
    vdg.getNumbers(length, numtri);
 
127
 
 
128
    dim[0] = length;
 
129
    dim[1] = 2;
 
130
    edge_db = PyArray_SimpleNew(2, dim, PyArray_INT);
 
131
    if (!edge_db) goto fail;
 
132
    edge_db_ptr = (int*)PyArray_DATA(edge_db);
 
133
    
 
134
    dim[0] = numtri;
 
135
    vertices = PyArray_SimpleNew(2, dim, PyArray_DOUBLE);
 
136
    if (!vertices) goto fail;
 
137
    vertices_ptr = (double*)PyArray_DATA(vertices);
 
138
 
 
139
    dim[1] = 3;
 
140
    tri_edges = PyArray_SimpleNew(2, dim, PyArray_INT);
 
141
    if (!tri_edges) goto fail;
 
142
    tri_edges_ptr = (int*)PyArray_DATA(tri_edges);
 
143
 
 
144
    tri_nbrs = PyArray_SimpleNew(2, dim, PyArray_INT);
 
145
    if (!tri_nbrs) goto fail;
 
146
    tri_nbrs_ptr = (int*)PyArray_DATA(tri_nbrs);
 
147
 
 
148
    for (i=0; i<(3*numtri); i++) {
 
149
        tri_edges_ptr[i] = tri_nbrs_ptr[i] = -1;
 
150
    }
 
151
 
 
152
    vdg.resetEdgeListIter();
 
153
    i = -1;
 
154
    while (vdg.getNextDelaunay(tri0, tri0x, tri0y, tri1, tri1x, tri1y, reg0, reg1)) {
 
155
        i++;
 
156
        INDEX2(edge_db_ptr,i,0) = reg0;
 
157
        INDEX2(edge_db_ptr,i,1) = reg1;
 
158
        if (tri0 > -1) {
 
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;
 
166
                    break;
 
167
                }
 
168
            }
 
169
        }
 
170
        if (tri1 > -1) {
 
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;
 
178
                    break;
 
179
                }
 
180
            }
 
181
        }
 
182
    }
 
183
 
 
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, 
 
188
        tri_nbrs_ptr);
 
189
 
 
190
    temp = Py_BuildValue("(OOOO)", vertices, edge_db, tri_edges, tri_nbrs);
 
191
    if (!temp) goto fail;
 
192
 
 
193
    Py_DECREF(vertices);
 
194
    Py_DECREF(edge_db);
 
195
    Py_DECREF(tri_edges);
 
196
    Py_DECREF(tri_nbrs);
 
197
 
 
198
    return temp;
 
199
 
 
200
fail:
 
201
    Py_XDECREF(vertices);
 
202
    Py_XDECREF(edge_db);
 
203
    Py_XDECREF(tri_edges);
 
204
    Py_XDECREF(tri_nbrs);
 
205
    return NULL;
 
206
}
 
207
 
 
208
static PyObject *linear_planes(int ntriangles, double *x, double *y, double *z,
 
209
    int *nodes)
 
210
{
 
211
    intp dims[2];
 
212
    PyObject *planes;
 
213
    int i;
 
214
    double *planes_ptr;
 
215
    double x02, y02, z02, x12, y12, z12, xy0212;
 
216
    
 
217
    dims[0] = ntriangles;
 
218
    dims[1] = 3;
 
219
    planes = PyArray_SimpleNew(2, dims, PyArray_DOUBLE);
 
220
    if (!planes) return NULL;
 
221
    planes_ptr = (double *)PyArray_DATA(planes);
 
222
 
 
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)];
 
230
 
 
231
        if (y12 != 0.0) {
 
232
            xy0212 = y02/y12;
 
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)]);
 
238
        } else {
 
239
            xy0212 = x02/x12;
 
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)]);
 
245
        }
 
246
    }
 
247
 
 
248
    return (PyObject*)planes;
 
249
}
 
250
 
 
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)
 
254
{
 
255
    double *planes_ptr;
 
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));
 
264
}
 
265
 
 
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)
 
270
{
 
271
    int ix, iy;
 
272
    double dx, dy, targetx, targety;
 
273
    int rowtri, coltri, tri;
 
274
    PyObject *z;
 
275
    double *z_ptr;
 
276
    intp dims[2];
 
277
 
 
278
    dims[0] = ysteps;
 
279
    dims[1] = xsteps;
 
280
    z = PyArray_SimpleNew(2, dims, PyArray_DOUBLE);
 
281
    if (!z) return NULL;
 
282
    z_ptr = (double*)PyArray_DATA(z);
 
283
 
 
284
    dx = (x1 - x0) / (xsteps-1);
 
285
    dy = (y1 - y0) / (ysteps-1);
 
286
 
 
287
    rowtri = 0;
 
288
    for (iy=0; iy<ysteps; iy++) {
 
289
        targety = y0 + dy*iy;
 
290
        rowtri = walking_triangles(rowtri, x0, targety, x, y, nodes, neighbors);
 
291
        tri = rowtri;
 
292
        for (ix=0; ix<xsteps; ix++) {
 
293
            targetx = x0 + dx*ix;
 
294
            INDEXN(z_ptr, xsteps, iy, ix) = linear_interpolate_single(
 
295
                targetx, targety,
 
296
                x, y, nodes, neighbors, planes, defvalue, tri, &coltri);
 
297
            if (coltri != -1) tri = coltri;
 
298
        }
 
299
    }
 
300
 
 
301
    return z;
 
302
}
 
303
 
 
304
static PyObject *compute_planes_method(PyObject *self, PyObject *args)
 
305
{
 
306
    PyObject *pyx, *pyy, *pyz, *pynodes;
 
307
    PyArrayObject *x, *y, *z, *nodes;
 
308
    int npoints, ntriangles;
 
309
 
 
310
    PyObject *planes;
 
311
 
 
312
    if (!PyArg_ParseTuple(args, "OOOO", &pyx, &pyy, &pyz, &pynodes)) {
 
313
        return NULL;
 
314
    }
 
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");
 
320
        goto fail;
 
321
    }
 
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");
 
325
        goto fail;
 
326
    }
 
327
    nodes = (PyArrayObject*)PyArray_ContiguousFromObject(pynodes, PyArray_INT, 2, 2);
 
328
    if (!nodes) {
 
329
        PyErr_SetString(PyExc_ValueError, "nodes must be a 2-D array of ints");
 
330
        goto fail;
 
331
    }
 
332
    ntriangles = PyArray_DIM(nodes, 0);
 
333
    if (PyArray_DIM(nodes, 1) != 3) {
 
334
        PyErr_SetString(PyExc_ValueError, "nodes must have shape (ntriangles, 3)");
 
335
        goto fail;
 
336
    }
 
337
 
 
338
    planes = linear_planes(ntriangles, (double*)PyArray_DATA(x), 
 
339
        (double*)PyArray_DATA(y), (double*)PyArray_DATA(z), (int*)PyArray_DATA(nodes));
 
340
 
 
341
    Py_DECREF(pyx);
 
342
    Py_DECREF(pyy);
 
343
    Py_DECREF(pyz);
 
344
    Py_DECREF(pynodes);
 
345
 
 
346
    return planes;
 
347
 
 
348
fail:
 
349
    Py_XDECREF(pyx);
 
350
    Py_XDECREF(pyy);
 
351
    Py_XDECREF(pyz);
 
352
    Py_XDECREF(pynodes);
 
353
    return NULL;
 
354
}
 
355
 
 
356
static PyObject *linear_interpolate_method(PyObject *self, PyObject *args)
 
357
{
 
358
    double x0, x1, y0, y1, defvalue;
 
359
    int xsteps, ysteps;
 
360
    PyObject *pyplanes, *pyx, *pyy, *pynodes, *pyneighbors, *grid;
 
361
    int npoints;
 
362
 
 
363
    PyArrayObject *planes, *x, *y, *nodes, *neighbors;
 
364
 
 
365
    if (!PyArg_ParseTuple(args, "ddiddidOOOOO", &x0, &x1, &xsteps, &y0, &y1, &ysteps,
 
366
           &defvalue, &pyplanes, &pyx, &pyy, &pynodes, &pyneighbors)) {
 
367
        return NULL;
 
368
    }
 
369
    x = (PyArrayObject*)PyArray_ContiguousFromObject(pyx, PyArray_DOUBLE, 1, 1);
 
370
    y = (PyArrayObject*)PyArray_ContiguousFromObject(pyy, PyArray_DOUBLE, 1, 1);
 
371
    if (!x || !y) {
 
372
        PyErr_SetString(PyExc_ValueError, "x,y must be 1-D arrays of floats");
 
373
        goto fail;
 
374
    }
 
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");
 
378
        goto fail;
 
379
    }
 
380
    planes = (PyArrayObject*)PyArray_ContiguousFromObject(pyplanes, PyArray_DOUBLE, 2, 2);
 
381
    if (!planes) {
 
382
        PyErr_SetString(PyExc_ValueError, "planes must be a 2-D array of floats");
 
383
        goto fail;
 
384
    }
 
385
    nodes = (PyArrayObject*)PyArray_ContiguousFromObject(pynodes, PyArray_INT, 2, 2);
 
386
    if (!nodes) {
 
387
        PyErr_SetString(PyExc_ValueError, "nodes must be a 2-D array of ints");
 
388
        goto fail;
 
389
    }
 
390
    neighbors = (PyArrayObject*)PyArray_ContiguousFromObject(pyneighbors, PyArray_INT, 2, 2);
 
391
    if (!neighbors) {
 
392
        PyErr_SetString(PyExc_ValueError, "neighbors must be a 2-D array of ints");
 
393
        goto fail;
 
394
    }
 
395
 
 
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));
 
400
    
 
401
    Py_DECREF(x);
 
402
    Py_DECREF(y);
 
403
    Py_DECREF(planes);
 
404
    Py_DECREF(nodes);
 
405
    Py_DECREF(neighbors);
 
406
 
 
407
    return grid;
 
408
 
 
409
fail:
 
410
    Py_XDECREF(x);
 
411
    Py_XDECREF(y);
 
412
    Py_XDECREF(planes);
 
413
    Py_XDECREF(nodes);
 
414
    Py_XDECREF(neighbors);
 
415
    return NULL;
 
416
}
 
417
 
 
418
#define CLEANUP \
 
419
    Py_XDECREF(x);\
 
420
    Py_XDECREF(y);\
 
421
    Py_XDECREF(z);\
 
422
    Py_XDECREF(intx);\
 
423
    Py_XDECREF(inty);\
 
424
    Py_XDECREF(centers);\
 
425
    Py_XDECREF(nodes);\
 
426
    Py_XDECREF(neighbors);\
 
427
    Py_XDECREF(intz);\
 
428
    return NULL;
 
429
 
 
430
static PyObject *nn_interpolate_unstructured_method(PyObject *self, PyObject *args)
 
431
{
 
432
    PyObject *pyx, *pyy, *pyz, *pycenters, *pynodes, *pyneighbors, *pyintx, *pyinty;
 
433
    PyArrayObject *x, *y, *z, *centers, *nodes, *neighbors, *intx, *inty, *intz;
 
434
    double defvalue;
 
435
    int size, npoints, ntriangles;
 
436
 
 
437
    if (!PyArg_ParseTuple(args, "OOdOOOOOO", &pyintx, &pyinty, &defvalue,
 
438
        &pyx, &pyy, &pyz, &pycenters, &pynodes, &pyneighbors)) {
 
439
        return NULL;
 
440
    }
 
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");
 
446
        //goto fail;
 
447
        CLEANUP
 
448
    }
 
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");
 
452
        //goto fail;
 
453
        CLEANUP
 
454
    }
 
455
    centers = (PyArrayObject*)PyArray_ContiguousFromObject(pycenters, PyArray_DOUBLE, 2, 2);
 
456
    if (!centers) {
 
457
        PyErr_SetString(PyExc_ValueError, "centers must be a 2-D array of ints");
 
458
        //goto fail;
 
459
        CLEANUP
 
460
    }
 
461
    nodes = (PyArrayObject*)PyArray_ContiguousFromObject(pynodes, PyArray_INT, 2, 2);
 
462
    if (!nodes) {
 
463
        PyErr_SetString(PyExc_ValueError, "nodes must be a 2-D array of ints");
 
464
        //goto fail;
 
465
        CLEANUP
 
466
    }
 
467
    neighbors = (PyArrayObject*)PyArray_ContiguousFromObject(pyneighbors, PyArray_INT, 2, 2);
 
468
    if (!neighbors) {
 
469
        PyErr_SetString(PyExc_ValueError, "neighbors must be a 2-D array of ints");
 
470
        //goto fail;
 
471
        CLEANUP
 
472
    }
 
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");
 
477
        //goto fail;
 
478
        CLEANUP
 
479
    }
 
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");
 
484
        CLEANUP
 
485
    }
 
486
    if (intx->nd != inty->nd) {
 
487
        PyErr_SetString(PyExc_ValueError, "intx,inty must have same shapes");
 
488
        CLEANUP
 
489
    }
 
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");
 
493
            CLEANUP
 
494
        }
 
495
    }
 
496
    intz = (PyArrayObject*)PyArray_SimpleNew(intx->nd, intx->dimensions, PyArray_DOUBLE);
 
497
    if (!intz) {CLEANUP}
 
498
 
 
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);
 
507
 
 
508
    return (PyObject*)intz;
 
509
}
 
510
 
 
511
#undef CLEANUP
 
512
 
 
513
#define CLEANUP \
 
514
    Py_XDECREF(x);\
 
515
    Py_XDECREF(y);\
 
516
    Py_XDECREF(z);\
 
517
    Py_XDECREF(centers);\
 
518
    Py_XDECREF(nodes);\
 
519
    Py_XDECREF(neighbors);\
 
520
    return NULL;
 
521
 
 
522
static PyObject *nn_interpolate_method(PyObject *self, PyObject *args)
 
523
{
 
524
    PyObject *pyx, *pyy, *pyz, *pycenters, *pynodes, *pyneighbors, *grid;
 
525
    PyArrayObject *x, *y, *z, *centers, *nodes, *neighbors;
 
526
    double x0, x1, y0, y1, defvalue;
 
527
    int xsteps, ysteps;
 
528
    int npoints, ntriangles;
 
529
    intp dims[2];
 
530
 
 
531
    if (!PyArg_ParseTuple(args, "ddiddidOOOOOO", &x0, &x1, &xsteps, 
 
532
        &y0, &y1, &ysteps, &defvalue, &pyx, &pyy, &pyz, &pycenters, &pynodes,
 
533
        &pyneighbors)) {
 
534
        return NULL;
 
535
    }
 
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");
 
541
        //goto fail;
 
542
        CLEANUP
 
543
    }
 
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");
 
547
        //goto fail;
 
548
        CLEANUP
 
549
    }
 
550
    centers = (PyArrayObject*)PyArray_ContiguousFromObject(pycenters, PyArray_DOUBLE, 2, 2);
 
551
    if (!centers) {
 
552
        PyErr_SetString(PyExc_ValueError, "centers must be a 2-D array of ints");
 
553
        //goto fail;
 
554
        CLEANUP
 
555
    }
 
556
    nodes = (PyArrayObject*)PyArray_ContiguousFromObject(pynodes, PyArray_INT, 2, 2);
 
557
    if (!nodes) {
 
558
        PyErr_SetString(PyExc_ValueError, "nodes must be a 2-D array of ints");
 
559
        //goto fail;
 
560
        CLEANUP
 
561
    }
 
562
    neighbors = (PyArrayObject*)PyArray_ContiguousFromObject(pyneighbors, PyArray_INT, 2, 2);
 
563
    if (!neighbors) {
 
564
        PyErr_SetString(PyExc_ValueError, "neighbors must be a 2-D array of ints");
 
565
        //goto fail;
 
566
        CLEANUP
 
567
    }
 
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");
 
572
        //goto fail;
 
573
        CLEANUP
 
574
    }
 
575
    //goto succeed; // XXX: Can't cross NaturalNeighbors instantiation with goto
 
576
 
 
577
    dims[0] = ysteps;
 
578
    dims[1] = xsteps;
 
579
    grid = PyArray_SimpleNew(2, dims, PyArray_DOUBLE);
 
580
    if (!grid) {CLEANUP} // goto fail;
 
581
 
 
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), 
 
587
        x0, x1, xsteps,
 
588
        y0, y1, ysteps,
 
589
        (double*)PyArray_DATA(grid),
 
590
        defvalue, 0);
 
591
 
 
592
    Py_DECREF(x);
 
593
    Py_DECREF(y);
 
594
    Py_DECREF(z);
 
595
    Py_DECREF(centers);
 
596
    Py_DECREF(nodes);
 
597
    Py_DECREF(neighbors);
 
598
 
 
599
    return grid;
 
600
 
 
601
}
 
602
#undef CLEANUP
 
603
 
 
604
static PyObject *delaunay_method(PyObject *self, PyObject *args)
 
605
{
 
606
    PyObject *pyx, *pyy, *mesh;
 
607
    PyArrayObject *x, *y;
 
608
    int npoints;
 
609
 
 
610
    if (!PyArg_ParseTuple(args, "OO", &pyx, &pyy)) {
 
611
        return NULL;
 
612
    }
 
613
    x = (PyArrayObject*)PyArray_ContiguousFromObject(pyx, PyArray_DOUBLE, 1, 1);
 
614
    y = (PyArrayObject*)PyArray_ContiguousFromObject(pyy, PyArray_DOUBLE, 1, 1);
 
615
    if (!x || !y) {
 
616
        PyErr_SetString(PyExc_ValueError, "x,y must be 1-D arrays of floats");
 
617
        goto fail;
 
618
    }
 
619
 
 
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");
 
623
        goto fail;
 
624
    }
 
625
 
 
626
    mesh = getMesh(npoints, (double*)(x->data), (double*)(y->data));
 
627
 
 
628
    if (!mesh) goto fail;
 
629
 
 
630
    Py_DECREF(x);
 
631
    Py_DECREF(y);
 
632
 
 
633
    return mesh;
 
634
 
 
635
fail:
 
636
    Py_XDECREF(x);
 
637
    Py_XDECREF(y);
 
638
    return NULL;
 
639
}
 
640
 
 
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,
 
655
        ""},
 
656
    {"linear_interpolate_grid", (PyCFunction)linear_interpolate_method, METH_VARARGS,
 
657
        ""},
 
658
    {"nn_interpolate_grid", (PyCFunction)nn_interpolate_method, METH_VARARGS,
 
659
        ""},
 
660
    {"nn_interpolate_unstructured", (PyCFunction)nn_interpolate_unstructured_method, METH_VARARGS,
 
661
        ""},
 
662
    {NULL, NULL, 0, NULL}
 
663
};
 
664
 
 
665
 
 
666
PyMODINIT_FUNC init_delaunay(void)
 
667
{
 
668
    PyObject* m;
 
669
    m = Py_InitModule3("_delaunay", delaunay_methods, 
 
670
        "Tools for computing the Delaunay triangulation and some operations on it.\n"
 
671
        );
 
672
    if (m == NULL)
 
673
        return;
 
674
    import_array();
 
675
}
 
676
 
 
677
} // extern "C"