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>
38 #ifndef SG_HAVE_NATIVE_SGI_COMPILERS
39 SG_USING_NAMESPACE(std);
44 static const double PI = 4.0*atan(1.0);
45 static const double TWOPI = 2.0*PI;
47 //---------------------------------------------------------------------------
48 SphereInterpolate::SphereInterpolate (int n, const double* x,
49 const double* y, const double* z,
50 const unsigned int* f)
52 // Assumes (x[i],y[i],z[i]) is unit length for all 0 <= i < n.
53 // For complete spherical coverage, include the two antipodal points
54 // (0,0,1,f(0,0,1)) and (0,0,-1,f(0,0,-1)) in the data set.
56 cout << "Initialising spherical interpolator.\n";
57 cout << "[ 0%] Allocating memory \r";
59 theta = new double[3*n];
60 phi = new double[3*n];
61 func = new unsigned int[3*n];
63 // convert data to spherical coordinates
66 for (i = 0; i < n; i++)
68 GetSphericalCoords(x[i],y[i],z[i],theta[i],phi[i]);
72 // use periodicity to get wrap-around in the Delaunay triangulation
73 cout << "[ 10%] copying vertices for wrap-around\r";
75 for (i = 0, j = n, k = 2*n; i < n; i++, j++, k++)
77 theta[j] = theta[i]+TWOPI;
78 theta[k] = theta[i]-TWOPI;
85 pInterp = new mgcLinInterp2D(3*n,theta,phi,func);
87 cout << "[100%] Finished initialising spherical interpolator. \n";
90 SphereInterpolate::SphereInterpolate (int n, const sgVec2* p, const unsigned int* f)
92 // Assumes (x[i],y[i],z[i]) is unit length for all 0 <= i < n.
93 // For complete spherical coverage, include the two antipodal points
94 // (0,0,1,f(0,0,1)) and (0,0,-1,f(0,0,-1)) in the data set.
95 cout << "Initialising spherical interpolator.\n";
96 cout << "[ 0%] Allocating memory \r";
98 theta = new double[3*n];
99 phi = new double[3*n];
100 func = new unsigned int[3*n];
102 // convert data to spherical coordinates
103 cout << "[ 10%] copying vertices for wrap-around \r";
106 for (i = 0, j = n, k = 2*n; i < n; i++, j++, k++)
112 // use periodicity to get wrap-around in the Delaunay triangulation
115 theta[j] = theta[i]+TWOPI;
116 theta[k] = theta[i]-TWOPI;
121 pInterp = new mgcLinInterp2D(3*n,theta,phi,func);
123 cout << "[100%] Finished initialising spherical interpolator. \n";
125 //---------------------------------------------------------------------------
126 SphereInterpolate::~SphereInterpolate ()
133 //---------------------------------------------------------------------------
134 void SphereInterpolate::GetSphericalCoords (const double x, const double y, const double z,
136 double& phiAngle) const
138 // Assumes (x,y,z) is unit length. Returns -PI <= thetaAngle <= PI
139 // and 0 <= phiAngle <= PI.
145 thetaAngle = atan2(y,x);
160 //---------------------------------------------------------------------------
161 int SphereInterpolate::Evaluate (const double x, const double y, const double z, EvaluateData& f) const
163 // assumes (x,y,z) is unit length
165 double thetaAngle, phiAngle;
166 GetSphericalCoords(x,y,z,thetaAngle,phiAngle);
167 return pInterp->Evaluate(thetaAngle,phiAngle,f);
169 //---------------------------------------------------------------------------
170 int SphereInterpolate::Evaluate (const double thetaAngle, const double phiAngle, EvaluateData& f) const
172 return pInterp->Evaluate(thetaAngle,phiAngle,f);
174 //---------------------------------------------------------------------------