2 WARNING - Do not remove this header.
4 This code is a templated version of the 'magic-software' spherical
5 interpolation code by Dave Eberly. The original (un-hacked) code can be
6 obtained from here: http://www.magic-software.com/gr_appr.htm
7 This code is derived from linintp2.h/cpp and sphrintp.h/cpp.
9 Dave Eberly says that the conditions for use are:
11 * You may distribute the original source code to others at no charge.
13 * You may modify the original source code and distribute it to others at
14 no charge. The modified code must be documented to indicate that it is
15 not part of the original package.
17 * You may use this code for non-commercial purposes. You may also
18 incorporate this code into commercial packages. However, you may not
19 sell any of your source code which contains my original and/or modified
20 source code. In such a case, you need to factor out my code and freely
23 * The original code comes with absolutely no warranty and no guarantee is
24 made that the code is bug-free.
26 This does not seem incompatible with GPL - so this modified version
27 is hereby placed under GPL along with the rest of FlightGear.
32 #include <simgear/compiler.h>
40 #ifndef SG_HAVE_NATIVE_SGI_COMPILERS
41 SG_USING_NAMESPACE(std);
45 //---------------------------------------------------------------------------
46 mgcLinInterp2D::mgcLinInterp2D (int _numPoints, double* x, double* y,
49 if ( (numPoints = _numPoints) < 3 )
58 cout << "[ 20%] allocating memory \r";
60 point = new double*[numPoints];
61 tmppoint = new double*[numPoints+3];
62 f = new unsigned int[numPoints];
64 for (i = 0; i < numPoints; i++)
65 point[i] = new double[2];
66 for (i = 0; i < numPoints+3; i++)
67 tmppoint[i] = new double[2];
68 for (i = 0; i < numPoints; i++)
70 point[i][0] = tmppoint[i][0] = x[i];
71 point[i][1] = tmppoint[i][1] = y[i];
76 cout << "[ 30%] creating delaunay diagram \r";
80 //---------------------------------------------------------------------------
81 mgcLinInterp2D::~mgcLinInterp2D ()
90 for (i = 0; i < numPoints; i++)
96 for (i = 0; i < numPoints+3; i++)
105 //---------------------------------------------------------------------------
106 void mgcLinInterp2D::ComputeBarycenter (Vertex& v0, Vertex& v1, Vertex& v2,
107 Vertex& ver, double c[3])
109 double a0 = v0.x-v2.x;
110 double b0 = v0.y-v2.y;
111 double a1 = v1.x-v2.x;
112 double b1 = v1.y-v2.y;
113 double a2 = ver.x-v2.x;
114 double b2 = ver.y-v2.y;
116 double m00 = a0*a0+b0*b0;
117 double m01 = a0*a1+b0*b1;
118 double m11 = a1*a1+b1*b1;
119 double r0 = a2*a0+b2*b0;
120 double r1 = a2*a1+b2*b1;
121 double det = m00*m11-m01*m01;
123 c[0] = (m11*r0-m01*r1)/det;
124 c[1] = (m00*r1-m01*r0)/det;
127 //---------------------------------------------------------------------------
128 int mgcLinInterp2D::InTriangle (Vertex& v0, Vertex& v1, Vertex& v2,
131 const double eps = 1e-08;
132 double tx, ty, nx, ny;
134 // test against normal to first edge
139 if ( tx*nx + ty*ny < -eps )
142 // test against normal to second edge
147 if ( tx*nx + ty*ny < -eps )
150 // test against normal to third edge
155 if ( tx*nx + ty*ny < -eps )
160 //---------------------------------------------------------------------------
161 int mgcLinInterp2D::Evaluate (double x, double y, EvaluateData& F)
163 Vertex ver = { x, y };
164 // determine which triangle contains the target point
168 for (i = 0; i < numTriangles; i++)
170 Triangle& t = triangle[i];
171 v0.x = point[t.vertex[0]][0];
172 v0.y = point[t.vertex[0]][1];
173 v1.x = point[t.vertex[1]][0];
174 v1.y = point[t.vertex[1]][1];
175 v2.x = point[t.vertex[2]][0];
176 v2.y = point[t.vertex[2]][1];
178 if ( InTriangle(v0,v1,v2,ver) )
182 if ( i == numTriangles ) // point is outside interpolation region
187 Triangle& t = triangle[i]; // (x,y) is in this triangle
189 // compute barycentric coordinates with respect to subtriangle
191 ComputeBarycenter(v0,v1,v2,ver,bary);
193 // compute barycentric combination of function values at vertices
194 F.index[0] = f[t.vertex[0]];
195 F.index[1] = f[t.vertex[1]];
196 F.index[2] = f[t.vertex[2]];
197 F.percentage[0] = bary[0];
198 F.percentage[1] = bary[1];
199 F.percentage[2] = bary[2];
203 //---------------------------------------------------------------------------
204 int mgcLinInterp2D::Delaunay2D ()
208 const double EPSILON = 1e-12;
209 const int TSIZE = 75;
210 const double RANGE = 10.0;
212 xmin = tmppoint[0][0];
214 ymin = tmppoint[0][1];
218 for (i = 0; i < numPoints; i++)
220 double value = tmppoint[i][0];
226 value = tmppoint[i][1];
233 double xrange = xmax-xmin, yrange = ymax-ymin;
234 double maxrange = xrange;
235 if ( maxrange < yrange )
238 // need to scale the data later to do a correct triangle count
239 double maxrange2 = maxrange*maxrange;
241 // tweak the points by very small random numbers
242 double bgs = EPSILON*maxrange;
244 for (i = 0; i < numPoints; i++)
246 tmppoint[i][0] += bgs*(0.5 - rand()/double(RAND_MAX));
247 tmppoint[i][1] += bgs*(0.5 - rand()/double(RAND_MAX));
252 { 5*RANGE, -RANGE, -RANGE },
253 { -RANGE, 5*RANGE, -RANGE }
255 for (i = 0; i < 3; i++)
257 tmppoint[numPoints+i][0] = xmin+xrange*wrk[0][i];
258 tmppoint[numPoints+i][1] = ymin+yrange*wrk[1][i];
261 int i0, i1, i2, i3, i4, i5, i6, i7, i8, i9, i11;
266 int** tmp = new int*[tsz+1];
267 tmp[0] = new int[2*(tsz+1)];
268 for (i0 = 1; i0 < tsz+1; i0++)
269 tmp[i0] = tmp[0] + 2*i0;
270 i1 = 2*(numPoints + 2);
272 int* id = new int[i1];
273 for (i0 = 0; i0 < i1; i0++)
276 int** a3s = new int*[i1];
277 a3s[0] = new int[3*i1];
278 for (i0 = 1; i0 < i1; i0++)
279 a3s[i0] = a3s[0] + 3*i0;
280 a3s[0][0] = numPoints;
281 a3s[0][1] = numPoints+1;
282 a3s[0][2] = numPoints+2;
284 double** ccr = new double*[i1]; // circumscribed centers and radii
285 ccr[0] = new double[3*i1];
286 for (i0 = 1; i0 < i1; i0++)
287 ccr[i0] = ccr[0] + 3*i0;
292 nts = 1; // number of triangles
295 cout << "[ 40%] create triangulation \r";
297 // compute triangulation
298 for (i0 = 0; i0 < numPoints; i0++)
302 for (i11 = 0; i11 < nts; i11++)
305 while ( a3s[i1][0] < 0 )
308 for (i2 = 0; i2 < 2; i2++)
310 double z = tmppoint[i0][i2]-ccr[i1][i2];
318 for (i2 = 0; i2 < 3; i2++)
323 for (i3 = 1; i3 < 2; i3++)
325 ii[i3] = ii[i3-1] + 1;
332 for (i3 = 0; i3 <= i8; i3++)
334 for (i5 = 0; i5 < 2; i5++)
335 if ( a3s[i1][ii[i5]] != tmp[i3][i5] )
337 for (i6 = 0; i6 < 2; i6++)
338 tmp[i3][i6] = tmp[i8][i6];
346 // temporary storage exceeded, increase TSIZE
350 for (i3 = 0; i3 < 2; i3++)
351 tmp[i7][i3] = a3s[i1][ii[i3]];
358 for (i1 = 0; i1 <= i7; i1++)
360 for (i2 = 0; i2 < 2; i2++)
361 for (wrk[i2][2] = 0, i3 = 0; i3 < 2; i3++)
363 wrk[i2][i3] = tmppoint[tmp[i1][i2]][i3]-tmppoint[i0][i3];
365 0.5*wrk[i2][i3]*(tmppoint[tmp[i1][i2]][i3]+
369 xx = wrk[0][0]*wrk[1][1]-wrk[1][0]*wrk[0][1];
370 ccr[id[i4]][0] = (wrk[0][2]*wrk[1][1]-wrk[1][2]*wrk[0][1])/xx;
371 ccr[id[i4]][1] = (wrk[0][0]*wrk[1][2]-wrk[1][0]*wrk[0][2])/xx;
373 for (ccr[id[i4]][2] = 0, i2 = 0; i2 < 2; i2++)
375 double z = tmppoint[i0][i2]-ccr[id[i4]][i2];
376 ccr[id[i4]][2] += z*z;
377 a3s[id[i4]][i2] = tmp[i1][i2];
387 // count the number of triangles
388 cout << "[ 50%] count the number of triangles \r";
392 for (i11 = 0; i11 < nts; i11++)
395 while ( a3s[i0][0] < 0 )
397 if ( a3s[i0][0] < numPoints )
399 for (i1 = 0; i1 < 2; i1++)
400 for (i2 = 0; i2 < 2; i2++)
402 tmppoint[a3s[i0][i1]][i2]-tmppoint[a3s[i0][2]][i2];
404 if ( fabs(wrk[0][0]*wrk[1][1]-wrk[0][1]*wrk[1][0]) > EPSILON*maxrange2 )
409 // create the triangles
410 cout << "[ 60%] create the triangles \r";
412 triangle = new Triangle[numTriangles];
416 for (i11 = 0; i11 < nts; i11++)
419 while ( a3s[i0][0] < 0 )
421 if ( a3s[i0][0] < numPoints )
423 for (i1 = 0; i1 < 2; i1++)
424 for (i2 = 0; i2 < 2; i2++)
426 tmppoint[a3s[i0][i1]][i2]-tmppoint[a3s[i0][2]][i2];
427 xx = wrk[0][0]*wrk[1][1]-wrk[0][1]*wrk[1][0];
428 if ( fabs(xx) > EPSILON*maxrange2 )
430 int delta = xx < 0 ? 1 : 0;
431 Triangle& tri = triangle[numTriangles];
432 tri.vertex[0] = a3s[i0][0];
433 tri.vertex[1] = a3s[i0][1+delta];
434 tri.vertex[2] = a3s[i0][2-delta];
444 cout << "[ 70%] build the edge table \r";
447 edge = new Edge[3*numTriangles];
450 for (i = 0; i < numTriangles; i++)
453 cout << "[ 7" << 10*i/numTriangles << "%] build the edge table \r";
455 Triangle& t = triangle[i];
457 for (j0 = 0, j1 = 1; j0 < 3; j0++, j1 = (j1+1)%3)
459 for (j = 0; j < numEdges; j++)
462 if ( (t.vertex[j0] == e.vertex[0]
463 && t.vertex[j1] == e.vertex[1])
464 || (t.vertex[j0] == e.vertex[1]
465 && t.vertex[j1] == e.vertex[0]) )
468 if ( j == numEdges ) // add edge to table
470 edge[j].vertex[0] = t.vertex[j0];
471 edge[j].vertex[1] = t.vertex[j1];
472 edge[j].triangle[0] = i;
473 edge[j].index[0] = j0;
474 edge[j].triangle[1] = -1;
477 else // edge already exists, add triangle to table
479 edge[j].triangle[1] = i;
480 edge[j].index[1] = j0;
485 // establish links between adjacent triangles
486 cout << "[ 80%] establishing links between adjacent triangles \r";
488 for (i = 0; i < numEdges; i++)
490 if ( edge[i].triangle[1] != -1 )
492 j0 = edge[i].triangle[0];
493 j1 = edge[i].triangle[1];
494 triangle[j0].adj[edge[i].index[0]] = j1;
495 triangle[j1].adj[edge[i].index[1]] = j0;
510 cout << "[ 90%] finsishes delauney triangulation \r";
514 //---------------------------------------------------------------------------
515 void mgcLinInterp2D::GetPoint (int i, double& x, double& y)
517 // assumes i is valid [can use PointCount() before passing i]
521 //---------------------------------------------------------------------------
522 void mgcLinInterp2D::GetEdge (int i, double& x0, double& y0, double& x1,
525 // assumes i is valid [can use EdgeCount() before passing i]
526 int v0 = edge[i].vertex[0], v1 = edge[i].vertex[1];
533 //---------------------------------------------------------------------------
534 void mgcLinInterp2D::GetTriangle (int i, double& x0, double& y0, double& x1,
535 double& y1, double& x2, double& y2)
537 // assumes i is valid [can use TriangleCount() before passing i]
538 int v0 = triangle[i].vertex[0];
539 int v1 = triangle[i].vertex[1];
540 int v2 = triangle[i].vertex[2];
549 //---------------------------------------------------------------------------