三次参数样条曲线拟合（任意控制点）
三次参数样条曲线拟合，主要是为解决三次样条函数不能解决的问题⽽提出的。三次样条函数要求x满⾜单调递增，即x序列满⾜
x0<x1<x2<...<xn。物理上的意义是，曲线不可以出现绕回或打圈。三次参数样条采⽤参数⽅程表⽰曲线，较为⽅便解决此问题。
三次样条函数原理百度⾥⽐较多，这⾥就不讲了。三次参数样条实现的原理是在三次样条函数的基础上进⾏参数化求解的，采⽤弦长累加⽅法。假设平⾯上控制点坐标P（xi，yi）i=0,1,2,...n。各点处累
加弦长可表⽰为：
.
因此，点P（x，y）坐标分量可表⽰为参数形式 x=x(s) 和 y=y(s) 。得到⼀张以累加弦长s作为参数的数据表：
这样就可以分别以s为变量，x和y为函数利⽤三次样条函数进⾏拟合。
三次样条函数实现代码：
#pragma once
class CSpline
{
public:
CSpline(void);
~CSpline(void);
void Splint(double *xa,double *ya,double *m,int n,double &x,double &y);
void Spline1(double *xa,double *ya,int n,double *&m,double bound1,double bound2);
void Spline2(double *xa,double *ya,int n,double *&m,double bound1=0,double bound2=0);
};
#include "StdAfx.h"
#include "3OrderSpline.h"
#include "math.h"
CSpline::CSpline(void)
{
}
CSpline::~CSpline(void)
{
}
//================================================================
// 函数功能：利⽤求出的⼆阶导数求给定点值(结合Spline1,Spline2)
// 输⼊参数： *xa 为横坐标值，ya为纵坐标值，n为点个数，m为⼆阶偏导数
/
/ x为给定点，y接收插出来的值
// 返回值：⽆返回值
//
// 作者：蒋锦朋 1034378054@qq
// 单位：中国地质⼤学（武汉）地球物理与空间信息学院
// ⽇期： 2014/12/03
//================================================================
void CSpline::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=abs(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为边界点⼀阶偏导数
// 返回值：⽆返回值
//
// 作者：蒋锦朋 1034378054@qq
// 单位：中国地质⼤学（武汉）地球物理与空间信息学院
// ⽇期： 2014/12/03
//=========================================================================== void CSpline::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];
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
// 单位：中国地质⼤学（武汉）地球物理与空间信息学院
// ⽇期： 2014/12/03
//=========================================================================== void CSpline::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;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;
}
#include "stdafx.h"
#include "stdio.h"
#include "math.h"
#include "3OrderSpline.h"
int _tmain(int argc, _TCHAR* argv[])
{
int num=8;
// initialize data
double x[8]={9.59,60.81,105.57,161.59,120.5,100.1,50.0,10.0};
double y[8]={61.97,107.13,56.56,105.27,120.5,150.0,110.0,180.0};
// generate parameterized variable
double *s=new double[num];
double *s=new double[num];
s[0]=0;
double tmp_s = 0;
for(int i=0;i<num-1;i++){
tmp_s=tmp_s+sqrt((x[i+1]-x[i])*(x[i+1]-x[i])+(y[i+1]-y[i])*(y[i+1]-y[i]));
s[i+1]=tmp_s;
}
// smoothed spline points
int N=tmp_s;
double *sd=new double[N];
double *xd=new double[N];
double *yd=new double[N];
for(int i=0;i<N-1;i++)
sd[i]=i;
sd[N-1]=tmp_s;
// ⾃然边界条件
//double f11x=0,f11y=0,f22x=0,f22y=0;
//CSpline spline;
//double *mx,*my;
//spline.Spline2(s,x,num,mx,f11x,f22x);
//spline.Spline2(s,y,num,my,f11y,f22y);
⼀阶导数边界条件（I型边界）
//double f1x=1,f1y=0,f2x=0,f2y=1; // ⾸端⼀阶导数(f1x,f1y)=(1,0)，尾端⼀阶导数(f1x,f1y)=(0,1)，其中(0,1)表⽰导数值⽆穷⼤
//CSpline spline;
//double *mx,*my;
//spline.Spline1(s,x,num,mx,f1x,f2x);s parameter
//spline.Spline1(s,y,num,my,f1y,f2y);
⼆阶导数边界条件（II型边界）
double f11x=-0.707,f11y=0.707,f22x=0,f22y=1; // ⾸端⼆阶导数(f11x,f11y)=(-0.707,0.707)，尾端⼆阶导数(f11x,f11y)=(0,1)，其中(0,1)表⽰导数值⽆穷⼤ CSpline spline;
double *mx,*my;
spline.Spline2(s,x,num,mx,f11x,f22x);
spline.Spline2(s,y,num,my,f11y,f22y);
for(int i=0;i<N;i++)
{
spline.Splint(s,x,mx,num,sd[i],xd[i]);
spline.Splint(s,y,my,num,sd[i],yd[i]);
}
FILE *fp_m_x = fopen("spline_", "wt");
FILE *fp_m_y = fopen("spline_", "wt");
for (int i = 0; i < N; i++){
fprintf(fp_m_x, "%lf\n", xd[i]);
fprintf(fp_m_y, "%lf\n", yd[i]);
}
fclose(fp_m_x);
fclose(fp_m_y);
delete sd;
delete xd;
delete yd;
delete mx;
delete my;
return 0;
}
程序运⾏后会⽣成拟合曲线坐标，⽣成的⽂件spline_和spline_分别为点的x和y坐标。利⽤matlab绘制得到如下图像：
版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。
发表评论