三次样条插值函数c语言（c++三次样条插值）

自己是一名石油工人，最近在工作过程中需要对一条处理出来的离散点测曲线进行插值，开始用的线性插值，插值出来的曲线在这些点处并不光滑（数学上叫一阶导数不连续），于是想到了用三次样条

最终在Numerical recipes in c这本书中找到了实现方法，试了一下，效果不错，其中前两个函数是三次样条的实现算法，第三个函数是调用方法，我都做了相应的注释，希望可以帮助大家解决自己工作学习过程中碰到的问题。

//求单点处的二阶导数，x、y、n为已知离散点横纵坐标及个数，yp1,ypn为一阶、二阶导数初始化边界条件，一般可以取0，y2为离散点处输出的二阶导数。
void spline(double *x, double *y, int n, double yp1, double ypn, double *y2)
{
	int i, k;
	--x;
	--y;
	--y2;
	double p, qn, sig, un, *u;
	u = new double[n + 1];
	if (yp1>0.99e30) y2[1] = u[1] = 0.0;
	else
	{
		y2[1] = -0.5;
		u[1] = (3.0 / (x[2] - x[1]))*((y[2] - y[1]) / (x[2] - x[1]) - yp1);
	}
	for (i = 2; i <= n - 1; i++)
	{
		sig = (x[i] - x[i - 1]) / (x[i + 1] - x[i - 1]);
		p = sig*y2[i - 1] + 2.0;
		y2[i] = (sig - 1.0) / p;
		u[i] = (y[i + 1] - y[i]) / (x[i + 1] - x[i]) - (y[i] - y[i - 1]) / (x[i] - x[i - 1]);
		u[i] = (6.0*u[i] / (x[i + 1] - x[i - 1]) - sig*u[i - 1]) / p;
	}
	if (ypn>0.99e30) qn = un = 0.0;
	else
	{
		qn = 0.5;
		un = (3.0 / (x[n] - x[n - 1]))*(ypn - (y[n] - y[n - 1]) / (x[n] - x[n - 1]));
	}
	y2[n] = (un - qn*u[n - 1]) / (qn*y2[n - 1] + 1.0);
	for (k = n - 1; k >= 1; k--)
		y2[k] = y2[k] * y2[k + 1] + u[k];
	delete u;
}

//计算任意点xx处的三次样条插值结果，x、y、n同上，y2a为上部函数计算的二阶导，xx为自变量。
double splint(double *x, double *y, int n, double *y2a, double xx)
{
	int k;
	double h, b, a;
	k = int((xx - x[0]) / (x[1] - x[0]));
	if (k>(n - 2) || k<0) return 0;
	h = x[k + 1] - x[k];
	if (h == 0.0) printf("Bad x input to routine splint");
	a = (x[k + 1] - xx) / h;
	b = (xx - x[k]) / h;
	return a*y[k] + b*y[k + 1] + ((a*a*a - a)*y2a[k] + (b*b*b - b)*y2a[k + 1])*(h*h) / 6.0;
}

//计算整条曲线的插值结果，用完别忘了删yout指针
float* SplineSplint(float deplevel, CArray<CURVE_INFO,CURVE_INFO> &flowarray, int &onum)
{
	int datanum=flowarray.GetCount();
	double* x=new double[datanum];
	double* y=new double[datanum];
	double* y2a=new double[datanum];
	for(int i=0;i<datanum;i++)
	{
		x[i]=flowarray.ElementAt(i).Depth;
		y[i]=flowarray.ElementAt(i).AverageData;
	}
	spline(x,y,datanum,0.0f,0.0f,y2a);
	onum=(flowarray.ElementAt(datanum-1).Depth-flowarray.ElementAt(0).Depth)/deplevel+1;
	float dep=flowarray.ElementAt(0).Depth;
	float* yout=new float[onum];
	for(int i=0;i<onum;i++)
	{
		dep = flowarray.ElementAt(0).Depth + i*deplevel;
		yout[i]=splint(x,y,datanum,y2a,dep);
	}
	delete[] x;
	delete[] y;
	delete[] y2a;
	return yout;
}
其中曲线信息的定义CURVE_INFO在头文件中，上面的函数只用到了结构中的Depth和AverageData，也就是实际离散点的横纵坐标，大家可以根据需要改成自己的数据：
typedef struct tagCURVE_INFO
{
	WIS_CURVE CurveInfo;
	float Depth;
	int PointNum;
	float AverageData;
}CURVE_INFO;