阿克瑪插值(2D)
#ifndef INTERPOLATION_HEADER<?xml:namespace prefix = o ns = "urn:schemas-microsoft-com:office:office" />
#define INTERPOLATION_HEADER
#include<cmath>
#include<set>
#include<vector>
#include<algorithm>
using namespace std;
#define DELTA 1e-6
#define VERY_CLOSE(x,y) fabs((x)-(y))<DELTA
template<typename T>
class Point2D
{
public:
Point2D(T _X=0,T _Y=0)
:X(_X),Y(_Y)
{
}
T X,Y;
};
bool MyFunc(const Point2D<double>& P1, const Point2D<double>& P2)
{
return P1.X+DELTA<P2.X;
}
class Point2D_X_Less_Delta
{
public:
template<typename T>
bool operator()(const
{
return P1.X+DELTA<P2.X;
}
};
template<typename T>
ostream & operator <<(ostream &os,const Point2D<T>& P)
{
return os<<P.X<<"/t"<<P.Y<<endl;
}
template<typename T>
istream & operator
{
return is>>P.X>>P.Y;
}
class InterpolationSmooth
{
public:
InterpolationSmooth(){}
template<typename InputIterator>
InterpolationSmooth(InputIterator iter1,InputIterator iter2):
Pts(iter1,iter2)
{
}
void Insert(const Point2D<double> &P)
{
Pts.insert(P);
}
template<typename InputIterator>
void Insert(InputIterator iter1,InputIterator iter2)
{
Pts.insert(iter1,iter2);
}
double GetValue(double X)
{
if (X<(*Pts.begin()).X)
{
throw exception("Iuput is too small!");
}
set<Point2D<double>,Point2D_X_Less_Delta >::iterator iter=lower_bound(Pts.begin(),Pts.end(),Point2D<double>(X,0),MyFunc);
if (iter==Pts.begin())
{
iter++;
}
int k=distance(Pts.begin(),iter);
k--;
double Xkp1=(*iter).X;
iter--;
double s0=(*iter).Y;
double s1=G(k);
double s2=(U(k)*3-G(k)*2-G(k+1))/(Xkp1-(*iter).X);
double s3=(G(k)+G(k+1)-U(k)*2.0)/((Xkp1-(*iter).X)*(Xkp1-(*iter).X));
s1*=X-(*iter).X;
s2*=(X-(*iter).X)*(X-(*iter).X);
s3*=(X-(*iter).X)*(X-(*iter).X)*(X-(*iter).X);
return s0+s1+s2+s3;
}
double GetDerivativeValue(double X)
{
if (X<(*Pts.begin()).X)
{
throw exception("Iuput is too small!");
}
set<Point2D<double>,Point2D_X_Less_Delta >::iterator iter=lower_bound(Pts.begin(),Pts.end(),Point2D<double>(X,0),MyFunc);
if (iter==Pts.begin())
{
iter++;
}
int k=distance(Pts.begin(),iter);
k--;
double Xkp1=(*iter).X;
iter--;
double s1=G(k);
double s2=(U(k)*3-G(k)*2-G(k+1))/(Xkp1-(*iter).X);
double s3=(G(k)+G(k+1)-U(k)*2.0)/((Xkp1-(*iter).X)*(Xkp1-(*iter).X));
s2*=(X-(*iter).X)*2.0;
s3*=(X-(*iter).X)*(X-(*iter).X)*3.0;
return s1+s2+s3;
}
set<Point2D<double>,Point2D_X_Less_Delta >::iterator begin()
{
return Pts.begin();
}
set<Point2D<double>,Point2D_X_Less_Delta >::iterator end()
{
return Pts.end();
}
size_t size()
{
return Pts.size();
}
private:
double U(int _Pos)
{
int Num=(int)(Pts.size());
if(Num<=2 ||_Pos>Num) throw exception("size is too small!");
if (_Pos>=Num-1)
{
return U(_Pos-1)*2.0-U(_Pos-2);
}
if (_Pos<0)
{
return U(_Pos+1)*2.0-U(_Pos+2);
}
set<Point2D<double>,Point2D_X_Less_Delta >::iterator iter1=Pts.begin();
set<Point2D<double>,Point2D_X_Less_Delta >::iterator iter2=Pts.begin();
advance(iter1,_Pos);
advance(iter2,_Pos+1);
return ((*iter2).Y-(*iter1).Y)/((*iter2).X-(*iter1).X);
}
double G(int _Pos)
{
double <?xml:namespace prefix = st1 ns = "urn:schemas-microsoft-com:office:smarttags" />Uk=U(_Pos),Ukp1=U(_Pos+1),Ukm1=U(_Pos-1),Ukm2=U(_Pos-2);
if ( VERY_CLOSE(Uk,Ukp1) && VERY_CLOSE(Ukm1,Ukm2) )
{
return (Uk+Ukm1)/2.0;
}
return (fabs(Ukp1-Uk)*Ukm1+fabs(Ukm1-Ukm2)*Uk)/(fabs(Ukp1-Uk)+fabs(Ukm1-Ukm2));
}
set<Point2D<double>,Point2D_X_Less_Delta > Pts;
};
#endif