三次参数样条曲线拟合（任意控制点）

最新推荐文章于 2024-05-20 13:29:22 发布

jiangjjp2812

最新推荐文章于 2024-05-20 13:29:22 发布

阅读量1.2w

点赞数 13

分类专栏： C++ 文章标签：三次参数样条曲线实现

本文链接：https://blog.csdn.net/jiangjjp2812/article/details/100140675

版权

C++ 专栏收录该内容

9 篇文章 6 订阅

订阅专栏

三次参数样条曲线拟合，主要是为解决三次样条函数不能解决的问题而提出的。三次样条函数要求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.com
// 单位：    中国地质大学（武汉） 地球物理与空间信息学院
// 日期：    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.com
// 单位：    中国地质大学（武汉） 地球物理与空间信息学院
// 日期：    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];

	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 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-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];
	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);
	//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_test_x.txt", "wt");
	FILE *fp_m_y = fopen("spline_test_y.txt", "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_test_x.txt和spline_test_y.txt分别为点的x和y坐标。利用matlab绘制得到如下图像：

附matlab检验代码：

clear all;
clc;

load spline_test_x.txt;    % 导入c++计算的结果
load spline_test_y.txt;

control_point_x=[9.59,60.81,105.57,161.59,120.5,100.1,50.0,10.0];
control_point_y=[61.97,107.13,56.56,105.27,120.5,150.0,110.0,180.0];

xy=zeros(8,2);
xy(:,1)=control_point_x;
xy(:,2)=control_point_y;
s=zeros(8,1);
s(1)=0;
tmp=0;
for i=1:7
    tmp=tmp+sqrt((xy(i+1,1)-xy(i,1))^2+(xy(i+1,2)-xy(i,2))^2);
    s(i+1)=tmp;
end

ss=s(1):1:s(8);

%cs3=csape(s',xy','variational');             %  自然边界,二阶导数值为0
%cs3=csape(s',xy','complete',[[1;0],[0;1]]);     %  一阶边界条件
cs3=csape(s',xy','second',[[-0.707;0.707],[0;1]]);     %  二阶边界条件
yy=ppval(cs3,ss);

figure;
h1=plot(control_point_x,control_point_y,'ok'); hold on;
h2=plot(spline_test_x,spline_test_y,'-r','LineWidth',2);
h3=plot(yy(1,:),yy(2,:),'--k','LineWidth',2);
legend([h1 h2 h3],'控制点','c++三次参数样条拟合','matlab三次参数样条拟合');

set(gca,'FontSize',15); set(gca,'FontWeight','bold'); set(gca,'LineWidth',2);
set(get(gca,'XLabel'),'Fontsize',15); set(get(gca,'YLabel'),'FontSize',15);

jiangjjp2812

关注

13
点赞
踩
84

收藏

觉得还不错? 一键收藏
9
评论
三次参数样条曲线拟合（任意控制点）

三次参数样条曲线拟合，主要是为解决三次样条函数不能解决的问题而提出的。三次样条函数要求x满足单调递增，即x序列满足x0<x1<x2<...<xn。物理上的意义是，曲线不可以出现绕回或打圈。三次参数样条采用参数方程表示曲线，较为方便解决此问题。三次样条函数原理百度里比较多，这里就不讲了。三次参数样条实现的原理是在三次样条函数的基础上进行参数化求解的...
复制链接

扫一扫