井眼轨迹的三次样条插值（vs + QT + coin3d）

最新推荐文章于 2022-02-07 17:51:50 发布

weixin_30369041

最新推荐文章于 2022-02-07 17:51:50 发布

阅读量557

点赞数

原文链接：http://www.cnblogs.com/ling123/p/5368575.html

版权

　　井眼轨迹数据的测量值是离散的，根据某些测斜公式，我们可以计算出离散的三维的井眼轨迹坐标，但是真实的井眼轨迹是一条平滑的曲线，这就需要我们对测斜数据进行插值，使井眼轨迹变得平滑，我暂时决定使用三次样条进行插值。（但是插值出来的点，并不是真实的测量值，并不能真实的反映经验轨迹的实际情况，仅供分析使用）

　　三次样条函数：（函数是在网上找到的，测试可用）

　　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);
    

};

View Code

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;
}

View Code

调用方法

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);
}

View Code

仅仅记录测方法的使用过程

　在使用过程中，要对横坐标进行从小到大的排序，还有去除重复点，否则结果出错

转载于:https://www.cnblogs.com/ling123/p/5368575.html

weixin_30369041

关注

0
点赞
踩
4

收藏

觉得还不错? 一键收藏
0
评论
井眼轨迹的三次样条插值（vs + QT + coin3d）

　　井眼轨迹数据的测量值是离散的，根据某些测斜公式，我们可以计算出离散的三维的井眼轨迹坐标，但是真实的井眼轨迹是一条平滑的曲线，这就需要我们对测斜数据进行插值，使井眼轨迹变得平滑，我暂时决定使用三次样条进行插值。（但是插值出来的点，并不是真实的测量值，并不能真实的反映经验轨迹的实际情况，仅供分析使用）　　三次样条函数：（函数是在网上找到的，测试可用）　　ThreeSpline.h...
复制链接

扫一扫