最近需要用到插值，但是總覺得線性插值得出來的太過硬了，所以想看一下三次樣條曲線怎么做。關於算法和程序實現的文章已經有很多了。這一篇文章寫下來主要的目的是為了
- 幫助自己理解，固化
- 已有的代碼不是在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還好用。