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还好用。