最近需要用到插值，但是总觉得线性插值得出来的太过硬了，所以想看一下三次样条曲线怎么做。关于算法和程序实现的文章已经有很多了。这一篇文章写下来主要的目的是为了
- 帮助自己理解，固化
- 已有的代码不是在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个了，可以求解了
怎么求解呢……有兴趣的话可以看看这篇：
这个方程组用矩阵形式写出来以后呢，可以发现是能够用TDMA算法算出来的，管他TDMA算法是什么呢……名字什么的并不重要
对于TDMA算法感兴趣的可以看看这篇：
如果觉得看三次曲线算法很烦的话，至少还是应该把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还好用。
