← Back to branch summary

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

« back to all changes in this revision

Viewing changes to Lib/sandbox/delaunay/natneighbors.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:
Lib/sandbox/delaunay/natneighbors.cpp
1
 
 
2
 
#include "delaunay_utils.h"
3
 
#include "natneighbors.h"
4
 
 
5
 
#include <stack>
6
 
#include <list>
7
 
#include <vector>
8
 
#include <set>
9
 
#include <math.h>
10
 
#include <iostream>
11
 
#include <iterator>
12
 
 
13
 
using namespace std;
14
 
 
15
 
NaturalNeighbors::NaturalNeighbors(int npoints, int ntriangles, double *x, double *y,
16
 
    double *centers, int *nodes, int *neighbors)
17
 
{
18
 
    this->npoints = npoints;
19
 
    this->ntriangles = ntriangles;
20
 
    this->x = x;
21
 
    this->y = y;
22
 
    this->centers = centers;
23
 
    this->nodes = nodes;
24
 
    this->neighbors = neighbors;
25
 
 
26
 
    this->radii2 = new double[ntriangles];
27
 
    for (int i=0; i<ntriangles; i++) {
28
 
        double x2 = x[INDEX3(nodes,i,0)] - INDEX2(centers,i,0);
29
 
        x2 = x2*x2;
30
 
        double y2 = y[INDEX3(nodes,i,0)] - INDEX2(centers,i,1);
31
 
        y2 = y2*y2;
32
 
        this->radii2[i] = x2 + y2;
33
 
    }
34
 
 
35
 
}
36
 
 
37
 
NaturalNeighbors::~NaturalNeighbors()
38
 
{
39
 
    delete[] this->radii2;
40
 
}
41
 
 
42
 
int NaturalNeighbors::find_containing_triangle(double targetx, double targety, int start_triangle)
43
 
{
44
 
    int final_triangle;
45
 
    final_triangle = walking_triangles(start_triangle, targetx, targety, 
46
 
        x, y, nodes, neighbors);
47
 
    return final_triangle;
48
 
}
49
 
 
50
 
double NaturalNeighbors::interpolate_one(double *z, double targetx, double targety,
51
 
    double defvalue, int &start_triangle)
52
 
{
53
 
    int t = find_containing_triangle(targetx, targety, start_triangle);
54
 
    if (t == -1) return defvalue;
55
 
 
56
 
    start_triangle = t;
57
 
    vector<int> circumtri;
58
 
 
59
 
    circumtri.push_back(t);
60
 
    stack<int> stackA;
61
 
    stack<int> stackB;
62
 
    int tnew, i;
63
 
 
64
 
    for (i=0; i<3; i++) {
65
 
        tnew = INDEX3(this->neighbors, t, i);
66
 
        if (tnew != -1) {
67
 
            stackA.push(tnew);
68
 
            stackB.push(t);
69
 
        }
70
 
    }
71
 
    while (!stackA.empty()) {
72
 
        tnew = stackA.top();
73
 
        stackA.pop();
74
 
        t = stackB.top();
75
 
        stackB.pop();
76
 
        double d2 = (SQ(targetx - INDEX2(this->centers,tnew,0))
77
 
                   + SQ(targety - INDEX2(this->centers,tnew,1)));
78
 
        if ((this->radii2[tnew]-d2) > TOLERANCE_EPS) {
79
 
            // tnew is a circumtriangle of the target
80
 
            circumtri.push_back(tnew);
81
 
            for (i=0; i<3; i++) {
82
 
                int ti = INDEX3(this->neighbors, tnew, i);
83
 
                if ((ti != -1) && (ti != t)) {
84
 
                    stackA.push(ti);
85
 
                    stackB.push(tnew);
86
 
                }
87
 
            }
88
 
        }
89
 
    }
90
 
 
91
 
    vector<int>::iterator it;
92
 
    double f = 0.0;
93
 
    double A = 0.0;
94
 
    double tA=0.0, yA=0.0, cA=0.0; // Kahan summation temps for A
95
 
    double tf=0.0, yf=0.0, cf=0.0; // Kahan summation temps for f
96
 
 
97
 
    vector<int> edge;
98
 
    bool onedge = false;
99
 
    bool onhull = false;
100
 
 
101
 
    for (it = circumtri.begin(); it != circumtri.end(); it++) {
102
 
        int t = *it;
103
 
        double vx = INDEX2(this->centers, t, 0);
104
 
        double vy = INDEX2(this->centers, t, 1);
105
 
        vector<double> c(6);
106
 
        for (int i=0; i<3; i++) {
107
 
            int j = EDGE0(i);
108
 
            int k = EDGE1(i);
109
 
 
110
 
            if (!circumcenter(
111
 
                    this->x[INDEX3(this->nodes, t, j)],
112
 
                    this->y[INDEX3(this->nodes, t, j)],
113
 
                    this->x[INDEX3(this->nodes, t, k)],
114
 
                    this->y[INDEX3(this->nodes, t, k)],
115
 
                    targetx, targety,
116
 
                    INDEX2(c, i, 0), INDEX2(c, i, 1))) {
117
 
 
118
 
                // bail out with the appropriate values if we're actually on a 
119
 
                // node
120
 
                if ((fabs(targetx - this->x[INDEX3(this->nodes, t, j)]) < TOLERANCE_EPS)
121
 
                 && (fabs(targety - this->y[INDEX3(this->nodes, t, j)]) < TOLERANCE_EPS)) {
122
 
                    return z[INDEX3(this->nodes, t, j)];
123
 
                } else if ((fabs(targetx - this->x[INDEX3(this->nodes, t, k)]) < TOLERANCE_EPS)
124
 
                        && (fabs(targety - this->y[INDEX3(this->nodes, t, k)]) < TOLERANCE_EPS)) {
125
 
                    return z[INDEX3(this->nodes, t, k)];
126
 
                } else if (!onedge) {
127
 
                    onedge = true;
128
 
                    edge.push_back(INDEX3(this->nodes, t, j));
129
 
                    edge.push_back(INDEX3(this->nodes, t, k));
130
 
                    onhull = (INDEX3(neighbors, t, i) == -1);
131
 
                }
132
 
            }
133
 
        }
134
 
        for (int i=0; i<3; i++) {
135
 
            int j = EDGE0(i);
136
 
            int k = EDGE1(i);
137
 
            int q = INDEX3(this->nodes, t, i);
138
 
            double ati = 0.0;
139
 
 
140
 
            if (!onedge || ((edge[0] != q) && edge[1] != q)) {
141
 
                ati = signed_area(vx, vy, 
142
 
                                  INDEX2(c, j, 0), INDEX2(c, j, 1),
143
 
                                  INDEX2(c, k, 0), INDEX2(c, k, 1));
144
 
 
145
 
 
146
 
                yA = ati - cA;
147
 
                tA = A + yA;
148
 
                cA = (tA - A) - yA;
149
 
                A = tA;
150
 
 
151
 
                yf = ati*z[q] - cf;
152
 
                tf = f + yf;
153
 
                cf = (tf - f) - yf;
154
 
                f = tf;
155
 
            }
156
 
        }
157
 
    }
158
 
 
159
 
    // If we're on an edge, then the scheme of adding up triangles as above
160
 
    // doesn't work so well. We'll take care of these two nodes here.
161
 
    if (onedge) {
162
 
 
163
 
        // If we're on the convex hull, then the other nodes don't actually 
164
 
        // contribute anything, just the nodes for the edge we're on. The 
165
 
        // Voronoi "polygons" are infinite in extent.
166
 
        if (onhull) {
167
 
            double a = (hypot(targetx-x[edge[0]], targety-y[edge[0]]) / 
168
 
                        hypot(x[edge[1]]-x[edge[0]], y[edge[1]]-y[edge[0]]));
169
 
            return (1-a) * z[edge[0]] + a*z[edge[1]];
170
 
        }
171
 
 
172
 
        set<int> T(circumtri.begin(), circumtri.end());
173
 
        vector<int> newedges0; // the two nodes that edge[0] still connect to
174
 
        vector<int> newedges1; // the two nodes that edge[1] still connect to
175
 
        set<int> alltri0; // all of the circumtriangle edge[0] participates in
176
 
        set<int> alltri1; // all of the circumtriangle edge[1] participates in
177
 
        for (it = circumtri.begin(); it != circumtri.end(); it++) {
178
 
            for (int i=0; i<3; i++) {
179
 
                int ti = INDEX3(this->neighbors, *it, i);
180
 
                int j = EDGE0(i);
181
 
                int k = EDGE1(i);
182
 
                int q0 = INDEX3(this->nodes, *it, j);
183
 
                int q1 = INDEX3(this->nodes, *it, k);
184
 
 
185
 
                if ((q0 == edge[0]) || (q1 == edge[0])) alltri0.insert(*it);
186
 
                if ((q0 == edge[1]) || (q1 == edge[1])) alltri1.insert(*it);
187
 
 
188
 
                if (!T.count(ti)) {
189
 
                    // neighbor is not in the set of circumtriangles
190
 
                    if (q0 == edge[0]) newedges0.push_back(q1);
191
 
                    if (q1 == edge[0]) newedges0.push_back(q0);
192
 
                    if (q0 == edge[1]) newedges1.push_back(q1);
193
 
                    if (q1 == edge[1]) newedges1.push_back(q0);
194
 
                }
195
 
            }
196
 
        }
197
 
 
198
 
        set<int>::iterator sit;
199
 
 
200
 
        double cx, cy;
201
 
        ConvexPolygon poly0;
202
 
        ConvexPolygon poly1;
203
 
 
204
 
        if (edge[1] != newedges0[0]) {
205
 
            circumcenter(this->x[edge[0]], this->y[edge[0]],
206
 
                         this->x[newedges0[0]], this->y[newedges0[0]],
207
 
                         targetx, targety,
208
 
                         cx, cy);
209
 
            poly0.push(cx, cy);
210
 
        }
211
 
        if (edge[1] != newedges0[1]) {
212
 
            circumcenter(this->x[edge[0]], this->y[edge[0]],
213
 
                         this->x[newedges0[1]], this->y[newedges0[1]],
214
 
                         targetx, targety,
215
 
                         cx, cy);
216
 
            poly0.push(cx, cy);
217
 
        }
218
 
 
219
 
        if (edge[0] != newedges1[0]) {
220
 
            circumcenter(this->x[edge[1]], this->y[edge[1]],
221
 
                         this->x[newedges1[0]], this->y[newedges1[0]],
222
 
                         targetx, targety,
223
 
                         cx, cy);
224
 
            poly1.push(cx, cy);
225
 
        }
226
 
        if (edge[0] != newedges1[1]) {
227
 
            circumcenter(this->x[edge[1]], this->y[edge[1]],
228
 
                         this->x[newedges1[1]], this->y[newedges1[1]],
229
 
                         targetx, targety,
230
 
                         cx, cy);
231
 
            poly1.push(cx, cy);
232
 
        }
233
 
 
234
 
        for (sit = alltri0.begin(); sit != alltri0.end(); sit++) {
235
 
            poly0.push(INDEX2(this->centers, *sit, 0),
236
 
                       INDEX2(this->centers, *sit, 1)); 
237
 
        }
238
 
        for (sit = alltri1.begin(); sit != alltri1.end(); sit++) {
239
 
            poly1.push(INDEX2(this->centers, *sit, 0),
240
 
                       INDEX2(this->centers, *sit, 1)); 
241
 
        }
242
 
 
243
 
        double a0 = poly0.area();
244
 
        double a1 = poly1.area();
245
 
 
246
 
 
247
 
        f += a0*z[edge[0]];
248
 
        A += a0;
249
 
        f += a1*z[edge[1]];
250
 
        A += a1;
251
 
 
252
 
        // Anticlimactic, isn't it?
253
 
    }
254
 
 
255
 
    f /= A;
256
 
    return f;
257
 
}
258
 
 
259
 
void NaturalNeighbors::interpolate_grid(double *z, 
260
 
    double x0, double x1, int xsteps,
261
 
    double y0, double y1, int ysteps,
262
 
    double *output,
263
 
    double defvalue, int start_triangle)
264
 
{
265
 
    int i, ix, iy, rowtri, coltri, tri;
266
 
    double dx, dy, targetx, targety;
267
 
 
268
 
    dx = (x1 - x0) / (xsteps-1);
269
 
    dy = (y1 - y0) / (ysteps-1);
270
 
 
271
 
    rowtri = 0;
272
 
    i = 0;
273
 
    for (iy=0; iy<ysteps; iy++) {
274
 
        targety = y0 + dy*iy;
275
 
        rowtri = find_containing_triangle(x0, targety, rowtri);
276
 
        tri = rowtri;
277
 
        for (ix=0; ix<xsteps; ix++) {
278
 
            targetx = x0 + dx*ix;
279
 
            coltri = tri;
280
 
            INDEXN(output, xsteps, iy, ix) = interpolate_one(z, targetx, targety,
281
 
                defvalue, coltri);
282
 
            if (coltri != -1) tri = coltri;
283
 
        }
284
 
    }
285
 
}
286
 
 
287
 
void NaturalNeighbors::interpolate_unstructured(double *z, int size, 
288
 
    double *intx, double *inty, double *output, double defvalue)
289
 
{
290
 
    int i, tri1, tri2;
291
 
 
292
 
    tri1 = 0;
293
 
    tri2 = 0;
294
 
    for (i=0; i<size; i++) {
295
 
        tri2 = tri1;
296
 
        output[i] = interpolate_one(z, intx[i], inty[i], defvalue, tri2);
297
 
        if (tri2 != -1) tri1 = tri2;
298
 
    }
299
 
}