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.
37 //---------------------------------------------------------------------------
39 mgcLinInterp2D<T>::mgcLinInterp2D (int _numPoints, double* x, double* y,
42 if ( (numPoints = _numPoints) < 3 )
51 cout << "[ 20%] allocating memory \r";
53 point = new double*[numPoints];
54 tmppoint = new double*[numPoints+3];
57 for (i = 0; i < numPoints; i++)
58 point[i] = new double[2];
59 for (i = 0; i < numPoints+3; i++)
60 tmppoint[i] = new double[2];
61 for (i = 0; i < numPoints; i++)
63 point[i][0] = tmppoint[i][0] = x[i];
64 point[i][1] = tmppoint[i][1] = y[i];
69 cout << "[ 30%] creating delaunay diagram \r";
73 //---------------------------------------------------------------------------
75 mgcLinInterp2D<T>::~mgcLinInterp2D ()
84 for (i = 0; i < numPoints; i++)
90 for (i = 0; i < numPoints+3; i++)
99 //---------------------------------------------------------------------------
101 void mgcLinInterp2D<T>::ComputeBarycenter (Vertex& v0, Vertex& v1, Vertex& v2,
102 Vertex& ver, double c[3])
104 double A0 = v0.x-v2.x;
105 double B0 = v0.y-v2.y;
106 double A1 = v1.x-v2.x;
107 double B1 = v1.y-v2.y;
108 double A2 = ver.x-v2.x;
109 double B2 = ver.y-v2.y;
111 double m00 = A0*A0+B0*B0;
112 double m01 = A0*A1+B0*B1;
113 double m11 = A1*A1+B1*B1;
114 double r0 = A2*A0+B2*B0;
115 double r1 = A2*A1+B2*B1;
116 double det = m00*m11-m01*m01;
118 c[0] = (m11*r0-m01*r1)/det;
119 c[1] = (m00*r1-m01*r0)/det;
122 //---------------------------------------------------------------------------
124 int mgcLinInterp2D<T>::InTriangle (Vertex& v0, Vertex& v1, Vertex& v2,
127 const double eps = 1e-08;
128 double tx, ty, nx, ny;
130 // test against normal to first edge
135 if ( tx*nx + ty*ny < -eps )
138 // test against normal to second edge
143 if ( tx*nx + ty*ny < -eps )
146 // test against normal to third edge
151 if ( tx*nx + ty*ny < -eps )
156 //---------------------------------------------------------------------------
158 int mgcLinInterp2D<T>::Evaluate (double x, double y, T& F)
160 Vertex ver = { x, y };
161 // determine which triangle contains the target point
165 for (i = 0; i < numTriangles; i++)
167 Triangle& t = triangle[i];
168 v0.x = point[t.vertex[0]][0];
169 v0.y = point[t.vertex[0]][1];
170 v1.x = point[t.vertex[1]][0];
171 v1.y = point[t.vertex[1]][1];
172 v2.x = point[t.vertex[2]][0];
173 v2.y = point[t.vertex[2]][1];
175 if ( InTriangle(v0,v1,v2,ver) )
179 if ( i == numTriangles ) // point is outside interpolation region
184 Triangle& t = triangle[i]; // (x,y) is in this triangle
186 // compute barycentric coordinates with respect to subtriangle
188 ComputeBarycenter(v0,v1,v2,ver,bary);
190 // compute barycentric combination of function values at vertices
191 F = bary[0]*f[t.vertex[0]]+bary[1]*f[t.vertex[1]]+bary[2]*f[t.vertex[2]];
195 //---------------------------------------------------------------------------
197 int mgcLinInterp2D<T>::Delaunay2D ()
201 const double EPSILON = 1e-12;
202 const int TSIZE = 75;
203 const double RANGE = 10.0;
205 xmin = tmppoint[0][0];
207 ymin = tmppoint[0][1];
211 for (i = 0; i < numPoints; i++)
213 double value = tmppoint[i][0];
219 value = tmppoint[i][1];
226 double xrange = xmax-xmin, yrange = ymax-ymin;
227 double maxrange = xrange;
228 if ( maxrange < yrange )
231 // need to scale the data later to do a correct triangle count
232 double maxrange2 = maxrange*maxrange;
234 // tweak the points by very small random numbers
235 double bgs = EPSILON*maxrange;
237 for (i = 0; i < numPoints; i++)
239 tmppoint[i][0] += bgs*(0.5 - rand()/double(RAND_MAX));
240 tmppoint[i][1] += bgs*(0.5 - rand()/double(RAND_MAX));
245 { 5*RANGE, -RANGE, -RANGE },
246 { -RANGE, 5*RANGE, -RANGE }
248 for (i = 0; i < 3; i++)
250 tmppoint[numPoints+i][0] = xmin+xrange*wrk[0][i];
251 tmppoint[numPoints+i][1] = ymin+yrange*wrk[1][i];
254 int i0, i1, i2, i3, i4, i5, i6, i7, i8, i9, i11;
259 int** tmp = new int*[tsz+1];
260 tmp[0] = new int[2*(tsz+1)];
261 for (i0 = 1; i0 < tsz+1; i0++)
262 tmp[i0] = tmp[0] + 2*i0;
263 i1 = 2*(numPoints + 2);
265 int* id = new int[i1];
266 for (i0 = 0; i0 < i1; i0++)
269 int** a3s = new int*[i1];
270 a3s[0] = new int[3*i1];
271 for (i0 = 1; i0 < i1; i0++)
272 a3s[i0] = a3s[0] + 3*i0;
273 a3s[0][0] = numPoints;
274 a3s[0][1] = numPoints+1;
275 a3s[0][2] = numPoints+2;
277 double** ccr = new double*[i1]; // circumscribed centers and radii
278 ccr[0] = new double[3*i1];
279 for (i0 = 1; i0 < i1; i0++)
280 ccr[i0] = ccr[0] + 3*i0;
285 nts = 1; // number of triangles
288 cout << "[ 40%] create triangulation \r";
290 // compute triangulation
291 for (i0 = 0; i0 < numPoints; i0++)
295 for (i11 = 0; i11 < nts; i11++)
298 while ( a3s[i1][0] < 0 )
301 for (i2 = 0; i2 < 2; i2++)
303 double z = tmppoint[i0][i2]-ccr[i1][i2];
311 for (i2 = 0; i2 < 3; i2++)
316 for (i3 = 1; i3 < 2; i3++)
318 ii[i3] = ii[i3-1] + 1;
325 for (i3 = 0; i3 <= i8; i3++)
327 for (i5 = 0; i5 < 2; i5++)
328 if ( a3s[i1][ii[i5]] != tmp[i3][i5] )
330 for (i6 = 0; i6 < 2; i6++)
331 tmp[i3][i6] = tmp[i8][i6];
339 // temporary storage exceeded, increase TSIZE
343 for (i3 = 0; i3 < 2; i3++)
344 tmp[i7][i3] = a3s[i1][ii[i3]];
351 for (i1 = 0; i1 <= i7; i1++)
353 for (i2 = 0; i2 < 2; i2++)
354 for (wrk[i2][2] = 0, i3 = 0; i3 < 2; i3++)
356 wrk[i2][i3] = tmppoint[tmp[i1][i2]][i3]-tmppoint[i0][i3];
358 0.5*wrk[i2][i3]*(tmppoint[tmp[i1][i2]][i3]+
362 xx = wrk[0][0]*wrk[1][1]-wrk[1][0]*wrk[0][1];
363 ccr[id[i4]][0] = (wrk[0][2]*wrk[1][1]-wrk[1][2]*wrk[0][1])/xx;
364 ccr[id[i4]][1] = (wrk[0][0]*wrk[1][2]-wrk[1][0]*wrk[0][2])/xx;
366 for (ccr[id[i4]][2] = 0, i2 = 0; i2 < 2; i2++)
368 double z = tmppoint[i0][i2]-ccr[id[i4]][i2];
369 ccr[id[i4]][2] += z*z;
370 a3s[id[i4]][i2] = tmp[i1][i2];
380 // count the number of triangles
381 cout << "[ 50%] count the number of triangles \r";
385 for (i11 = 0; i11 < nts; i11++)
388 while ( a3s[i0][0] < 0 )
390 if ( a3s[i0][0] < numPoints )
392 for (i1 = 0; i1 < 2; i1++)
393 for (i2 = 0; i2 < 2; i2++)
395 tmppoint[a3s[i0][i1]][i2]-tmppoint[a3s[i0][2]][i2];
397 if ( fabs(wrk[0][0]*wrk[1][1]-wrk[0][1]*wrk[1][0]) > EPSILON*maxrange2 )
402 // create the triangles
403 cout << "[ 60%] create the triangles \r";
405 triangle = new Triangle[numTriangles];
409 for (i11 = 0; i11 < nts; i11++)
412 while ( a3s[i0][0] < 0 )
414 if ( a3s[i0][0] < numPoints )
416 for (i1 = 0; i1 < 2; i1++)
417 for (i2 = 0; i2 < 2; i2++)
419 tmppoint[a3s[i0][i1]][i2]-tmppoint[a3s[i0][2]][i2];
420 xx = wrk[0][0]*wrk[1][1]-wrk[0][1]*wrk[1][0];
421 if ( fabs(xx) > EPSILON*maxrange2 )
423 int delta = xx < 0 ? 1 : 0;
424 Triangle& tri = triangle[numTriangles];
425 tri.vertex[0] = a3s[i0][0];
426 tri.vertex[1] = a3s[i0][1+delta];
427 tri.vertex[2] = a3s[i0][2-delta];
437 cout << "[ 70%] build the edge table \r";
440 edge = new Edge[3*numTriangles];
443 for (i = 0; i < numTriangles; i++)
446 cout << "[ 7" << 10*i/numTriangles << "%] build the edge table \r";
448 Triangle& t = triangle[i];
450 for (j0 = 0, j1 = 1; j0 < 3; j0++, j1 = (j1+1)%3)
452 for (j = 0; j < numEdges; j++)
455 if ( (t.vertex[j0] == e.vertex[0]
456 && t.vertex[j1] == e.vertex[1])
457 || (t.vertex[j0] == e.vertex[1]
458 && t.vertex[j1] == e.vertex[0]) )
461 if ( j == numEdges ) // add edge to table
463 edge[j].vertex[0] = t.vertex[j0];
464 edge[j].vertex[1] = t.vertex[j1];
465 edge[j].triangle[0] = i;
466 edge[j].index[0] = j0;
467 edge[j].triangle[1] = -1;
470 else // edge already exists, add triangle to table
472 edge[j].triangle[1] = i;
473 edge[j].index[1] = j0;
478 // establish links between adjacent triangles
479 cout << "[ 80%] establishing links between adjacent triangles \r";
481 for (i = 0; i < numEdges; i++)
483 if ( edge[i].triangle[1] != -1 )
485 j0 = edge[i].triangle[0];
486 j1 = edge[i].triangle[1];
487 triangle[j0].adj[edge[i].index[0]] = j1;
488 triangle[j1].adj[edge[i].index[1]] = j0;
503 cout << "[ 90%] finsishes delauney triangulation \r";
507 //---------------------------------------------------------------------------
509 void mgcLinInterp2D<T>::GetPoint (int i, double& x, double& y)
511 // assumes i is valid [can use PointCount() before passing i]
515 //---------------------------------------------------------------------------
517 void mgcLinInterp2D<T>::GetEdge (int i, double& x0, double& y0, double& x1,
520 // assumes i is valid [can use EdgeCount() before passing i]
521 int v0 = edge[i].vertex[0], v1 = edge[i].vertex[1];
528 //---------------------------------------------------------------------------
530 void mgcLinInterp2D<T>::GetTriangle (int i, double& x0, double& y0, double& x1,
531 double& y1, double& x2, double& y2)
533 // assumes i is valid [can use TriangleCount() before passing i]
534 int v0 = triangle[i].vertex[0];
535 int v1 = triangle[i].vertex[1];
536 int v2 = triangle[i].vertex[2];
545 //---------------------------------------------------------------------------