标签:
最近需要用到插值,但是总觉得线性插值得出来的太过硬了,所以想看一下三次样条曲线怎么做。关于算法和程序实现的文章已经有很多了。这一篇文章写下来主要的目的是为了
我的理解,分段三次样条曲线求解就是:
怎么求解呢……有兴趣的话可以看看这篇:
http://www.cnblogs.com/xpvincent/archive/2013/01/26/2878092.html
这个方程组用矩阵形式写出来以后呢,可以发现是能够用TDMA算法算出来的,管他TDMA算法是什么呢……名字什么的并不重要
对于TDMA算法感兴趣的可以看看这篇:
http://www.cnblogs.com/xpvincent/archive/2013/01/25/2877411.html
如果觉得看三次曲线算法很烦的话,至少还是应该把TDMA算法看一下,要不然接下来看这些代码也是眼花花的
目前只是做二维的插值,三维的话,估计是xy平面来一次,然后xz平面再来一次?
下面放代码,因为还只是试验一下,所以代码写的不怎么好看……
float[] x={1.4f,3f,4f,5.5f,6.5f,7.4f}; float[] y={2.3f,4.4f,5.6f,8f,10f,10.4f}; float[] a; float[] b; float[] c; float[] d; float[] h; float[] xExt; float[] yExt;
这一段只是把6个点的信息写进去,abcd是TDMA算法和三次曲线插值过程中反复会用到的,所以直接在类里面声明好了,h是三次曲线插值算法里面会用到的各个线段之间x分量上的长度。
abcd以及xy都是6的Length,h是5
xExt和yExt是存储插值扩增以后的点坐标的。
static float[] TDMA(float [] ta,float [] tb,float [] tc, float [] tx) { int n=tx.Length; tc[0]=tc[0]/tb[0]; tx[0]=tx[0]/tb[0]; for(int i=1;i<n;i++){ float m=1/(tb[i]-ta[i]*tc[i-1]); tc[i]=tc[i]*m; tx[i]=(tx[i]-ta[i]*tx[i-1])*m; } for(int i=n-2;i>0;i--) { tx[i]=tx[i]-tc[i]*tx[i+1]; } return tx; }
TDMA算法,程序基本上是从引用的文章中拿过来的。
这个,如果没看过TDMA算法是什么的话,我也真不知道怎么说。
这里面的ta,tb,tc分别对应左边矩阵的左下,中间右上三个斜线,tx对应等号右侧向量。
写成这样是为了方便运算,但是直接看起来有点反直觉。还是那句话,先去看算法吧
三次曲线各段系数的计算:
private void Cinterp(float[] x, float[] y) { int n=x.Length; float[] m=new float[n]; h=new float[n-1]; a=new float[n]; b=new float[n]; c=new float[n]; d=new float[n]; for(int i=0;i<n-1;i++) { h[i]=x[i+1]-x[i]; } a[0]=0; b[0]=1; c[0]=0; d[0]=0; a[n-1]=0; b[n-1]=1; c[n-1]=0; d[n-1]=0; for(int i=1;i<n-1;i++) { a[i]=h[i-1]; b[i]=2*(h[i-1]+h[i]); c[i]=h[i]; d[i]=6*((y[i+1]-y[i])/h[i]-(y[i]-y[i-1])/h[i-1]); } m=TDMA(a,b,c,d); for(int i=0;i<n-1;i++) { a[i]=y[i]; b[i]=(y[i+1]-y[i])/h[i]-h[i]*m[i]/2-h[i]*(m[i+1]-m[i])/6; c[i]=m[i]/2; d[i]=(m[i+1]-m[i])/(6*h[i]); } }
其实我写的好像也不怎么简略。
这里面,引用TDMA函数之前的abcd是代表着求解矩阵里面的abcd
引用TDMA函数之后的abcd则是三次函数里面的系数,abcd……
我只是懒得再去声明四个新的数组而已……
另外各段长度h后面还用得到,所以在类里面已经声明了,m的话外面用不到,就只在method里面声明了。
这里系数算好了以后就是插值了。我这里是给每两个点中间再插入两个新的点
先在x方向插好:
private float[] Xinterp(float[] xin,float[] hin) { float[] xExt=new float[(xin.Length-1)*3+1]; int i=0; for(;i<(xin.Length-1);i++) { xExt[i*3]=xin[i]; xExt[i*3+1]=xin[i]+hin[i]/3; xExt[i*3+2]=xin[i]+hin[i]/3*2; } xExt[i*3]=xin[i]; return xExt; }
然后根据x方向的点算出来y,既然每每个线段插2个点,那自然从起点开始每三个点用一条方程算啦
private float[] Yinterp(float[] xex) { float[] yout=new float[xex.Length]; for(int i=0;i<xex.Length;i++) { int seg=(int)i/3; float h=xex[i]-xex[seg*3]; yout[i]=a[seg]+b[seg]*h+c[seg]*h*h+d[seg]*h*h*h; } yout[xex.Length-1]=y[y.Length-1]; return yout; }
那么接下来在Scene view里面看看效果:
start():
void Start () { Cinterp(x,y); for(int i=0;i<(x.Length-1);i++){ // Debug.Log(a[i]+","+b[i]+","+c[i]+","+d[i]); } xExt=Xinterp(x,h); yExt=Yinterp(xExt); }
void OnDrawGizmos() { // for(int i=0;i<5;i++){ // Gizmos.DrawLine(new Vector2(x[i],y[i]),new Vector2(x[i+1],y[i+1])); // } for(int i=0;i<6;i++){ Gizmos.DrawSphere(new Vector2(x[i],y[i]),0.1f); } for(int i=0;i<xExt.Length-1;i++){ Gizmos.DrawLine(new Vector2(xExt[i],yExt[i]),new Vector2(xExt[i+1],yExt[i+1])); } }
顺便说一下,OnDrawGizmos真是调试神奇啊,比Debug.DrawLine还好用。
标签:
原文地址:http://www.cnblogs.com/SiumingLearning/p/CubicInterpolateInUnity3D.html