Unity中三次样条插值曲线的实现
最近需要用到插值，但是总觉得线性插值得出来的太过硬了，所以想看一下三次样条曲线怎么做。关于算法和程序实现的文章已经有很多了。这一篇文章写下来主要的目的是为了
- 帮助自己理解，固化
- 已有的代码不是在unity平台上实现的，所以代码相对繁杂，这里进一步做简化
我的理解，分段三次样条曲线求解就是：
- 已知：n个点，n-1个三次方程（a+bx+cx^2+dx^3)，而这些三次方程2一阶和二阶导数连续，这些三次方程当然在已知点上也是连续的
- 一阶二阶导数连续，就是在中间的连接点上（去掉头尾总共n-2个），前后两个方程导数相等
- 方程在已知点上连续就更好理解了，一条线贯穿全部点嘛
- 这样的话已知的条件就有点：n个，一阶导数连续n-2个，二阶导数连续n-2个，点连续n-2个，总共是4n-6个
- 未知：n-1个三次方程的系数abcd，这样子就有4n-4个未知数了
- 已知比未知少了两个条件，所以就会有各种边界设定，我们就用自然边界好了，也就是假设边界上没有力让曲线弯曲，大致上就是像ps里面曲线的端点不加限制方向的样子吧。有了边界条件设定以后呢，已知和未知都是4n-4个了，可以求解了
怎么求解呢……有兴趣的话可以看看这篇：
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还好用。
