井眼轨迹的三次样条插值 （vs + QT + coin3d）
井眼轨迹数据的测量值是离散的，根据某些测斜公式，我们可以计算出离散的三维的井眼轨迹坐标，但是真实的井眼轨迹是一条平滑的曲线，这就需要我们对测斜数据进行插值，使井眼轨迹变得平滑，我暂时决定使用三次样条进行插值。（但是插值出来的点，并不是真实的测量值，并不能真实的反映经验轨迹的实际情况，仅供分析使用）
三次样条函数：（函数是在网上找到的，测试可用）
ThreeSpline.h
#pragma once class ThreeSpline { public: ThreeSpline(void); ~ThreeSpline(void); // 利用求出的二阶导数求给定点值(结合Spline1,Spline2) void Splint(double *xa,double *ya,double *m,int n,double &x,double &y); //对一系列点求二阶偏导数，点横坐标单调递增（I型边界）(结合Spline)一阶偏导数 void Spline1(double *xa,double *ya,int n,double *&m,double bound1,double bound2); //对一系列点求二阶偏导数，点横坐标单调递增（II型边界）(结合Spline)二阶偏导数 void Spline2(double *xa,double *ya,int n,double *&m,double bound1=0,double bound2=0); };
ThreeSpline.cpp
#include "ThreeSpline.h" ThreeSpline::ThreeSpline(void) { } ThreeSpline::~ThreeSpline(void) { } //================================================================ // 函数功能： 利用求出的二阶导数求给定点值(结合Spline1,Spline2) // 输入参数： *xa 为横坐标值，ya为纵坐标值，n为点个数，m为二阶偏导数 // x为给定点，y接收插出来的值 // 返回值： 无返回值 //================================================================ void ThreeSpline::Splint(double *xa,double *ya,double *m,int n,double &x,double &y) { int klo,khi,k; klo=0; khi=n-1; double hh,bb,aa; while(khi-klo>1) // 二分法查找x所在区间段 { k=(khi+klo)>>1; if(xa[k]>x) khi=k; else klo=k; } hh=xa[khi]-xa[klo]; aa=(xa[khi]-x)/hh; bb=(x-xa[klo])/hh; y=aa*ya[klo]+bb*ya[khi]+((aa*aa*aa-aa)*m[klo]+(bb*bb*bb-bb)*m[khi])*hh*hh/6.0; } //=========================================================================== // 函数功能： 对一系列点求二阶偏导数，点横坐标单调递增（I型边界）(结合Spline) // 输入参数： *xa 为横坐标值，ya为纵坐标值，n为点个数，m为二阶偏导数(输出值） // bound1、bound2为边界点一阶偏导数 // 返回值： 无返回值 // //=========================================================================== void ThreeSpline::Spline1(double *xa,double *ya,int n,double *&m,double bound1,double bound2) { // 追赶法解方程求二阶偏导数 double f1=bound1,f2=bound2; double *a=new double[n]; // a:稀疏矩阵最下边一串数 double *b=new double[n]; // b:稀疏矩阵最中间一串数 double *c=new double[n]; // c:稀疏矩阵最上边一串数 double *d=new double[n]; double *f=new double[n]; double *bt=new double[n]; double *gm=new double[n]; double *h=new double[n]; m=new double[n]; for(int i=0;i<n;i++) b[i]=2; // 中间一串数为2 for(int i=0;i<n-1;i++) h[i]=xa[i+1]-xa[i]; // 各段步长 for(int i=1;i<n-1;i++) a[i]=h[i-1]/(h[i-1]+h[i]); a[n-1]=1; c[0]=1; for(int i=1;i<n-1;i++) c[i]=h[i]/(h[i-1]+h[i]); for(int i=0;i<n-1;i++) f[i]=(ya[i+1]-ya[i])/(xa[i+1]-xa[i]); d[0]=6*(f[0]-f1)/h[0]; d[n-1]=6*(f2-f[n-2])/h[n-2]; for(int i=1;i<n-1;i++) d[i]=6*(f[i]-f[i-1])/(h[i-1]+h[i]); bt[0]=c[0]/b[0]; // 追赶法求解方程 for(int i=1;i<n-1;i++) bt[i]=c[i]/(b[i]-a[i]*bt[i-1]); gm[0]=d[0]/b[0]; for(int i=1;i<=n-1;i++) gm[i]=(d[i]-a[i]*gm[i-1])/(b[i]-a[i]*bt[i-1]); m[n-1]=gm[n-1]; for(int i=n-2;i>=0;i--) m[i]=gm[i]-bt[i]*m[i+1]; delete a; delete b; delete c; delete d; delete gm; delete bt; delete f; delete h; } //=========================================================================== // 函数功能： 对一系列点求二阶偏导数，点横坐标单调递增（II型边界）(结合Spline) // 输入参数： *xa 为横坐标值，ya为纵坐标值，n为点个数，m为二阶偏导数(输出值） // bound1、bound2为边界点二阶偏导数，当bound1和bound2不给值时则使用 // 默认值0，即自然边界 // 返回值： 无返回值 // // 作者： 蒋锦朋 1034378054@qq.com // 单位： 中国地质大学（武汉） 地球物理与空间信息学院 // 日期： 2014/12/03 //=========================================================================== void ThreeSpline::Spline2(double *xa,double *ya,int n,double *&m,double bound1,double bound2) { // 追赶法解方程求二阶偏导数 double f11=bound1,f22=bound2; double *a=new double[n]; // a:稀疏矩阵最下边一串数 double *b=new double[n]; // b:稀疏矩阵最中间一串数 double *c=new double[n]; // c:稀疏矩阵最上边一串数 double *d=new double[n]; double *f=new double[n]; double *bt=new double[n]; double *gm=new double[n]; double *h=new double[n]; m=new double[n]; for(int i=0;i<n;i++) b[i]=2; for(int i=0;i<n-1;i++) h[i]=xa[i+1]-xa[i]; for(int i=1;i<n-1;i++) a[i]=h[i-1]/(h[i-1]+h[i]); a[n-1]=1; c[0]=1; for(int i=1;i<n-1;i++) c[i]=h[i]/(h[i-1]+h[i]); for(int i=0;i<n-1;i++) f[i]=(ya[i+1]-ya[i])/(xa[i+1]-xa[i]); //d[0]=6*(f[0]-f1)/h[0]; //d[n-1]=6*(f2-f[n-2])/h[n-2]; for(int i=1;i<n-1;i++) d[i]=6*(f[i]-f[i-1])/(h[i-1]+h[i]); d[1]=d[1]-a[1]*f11; d[n-2]=d[n-2]-c[n-2]*f22; //bt[0]=c[0]/b[0]; //for(int i=1;i<n-1;i++) bt[i]=c[i]/(b[i]-a[i]*bt[i-1]); //gm[0]=d[0]/b[0]; //for(int i=1;i<=n-1;i++) gm[i]=(d[i]-a[i]*gm[i-1])/(b[i]-a[i]*bt[i-1]); //m[n-1]=gm[n-1]; //for(int i=n-2;i>=0;i--) m[i]=gm[i]-bt[i]*m[i+1]; // 追赶法求解方程 bt[1]=c[1]/b[1]; for(int i=2;i<n-2;i++) bt[i]=c[i]/(b[i]-a[i]*bt[i-1]); gm[1]=d[1]/b[1]; for(int i=2;i<=n-2;i++) gm[i]=(d[i]-a[i]*gm[i-1])/(b[i]-a[i]*bt[i-1]); m[n-2]=gm[n-2]; for(int i=n-3;i>=1;i--) m[i]=gm[i]-bt[i]*m[i+1]; m[0]=f11; m[n-1]=f22; delete a; delete b; delete c; delete d; delete gm; delete bt; delete f; delete h; }
调用方法
void CMathTestView::OnSpline() { // TODO: 在此添加命令处理程序代码 double x[12]={0.52,8,17.95,28.65,50.65,104.6,156.6,260.7,364.4,468,507,520}; double y[12]={5.28794,13.84,20.2,24.9,31.1,36.5,36.6,31,20.9,7.8,1.5,0.2}; double xd[8]={4,14,30,60,130,230,450,515}; double f1=1.86548,f2=-0.046115; double f11=-0.279319,f22=0.0111560; CSpline spline; double *m; spline.Spline2(x,y,12,m,f11,f22); CString showstr; showstr.Empty(); for(int i=0;i<8;i++) { double yd; spline.Splint(x,y,m,12,xd[i],yd); CString s; s.Format(_T("%lf %lf\n"),xd[i],yd); showstr+=s; } delete m; MessageBox(showstr); }
仅仅记录测方法的使用过程
在使用过程中，要对横坐标进行从小到大的排序，还有去除重复点，否则结果出错
