样条之Akima光滑插值函数

核心代码：

 //////////////////////////////////////////////////////////////////////

 // Akima光滑插值

 // t    - 存放指定的插值点的值

 // s[]  - 一维数组，长度为5，其中s(0)，s(1)，s(2)，s(3)返回三次多项式的系数，

 //        s(4)返回指定插值点t处的函数近似值f(t)（k<0时）或任意值（k>=0时）

 // k    - 控制参数，若k>=0，则只计算第k个子区间[x(k), x(k+1)]上的三次多项式的系数

 //////////////////////////////////////////////////////////////////////

 static float GetValueAkima(const void* valuesPtr, int stride, int n, float t, float* s, int k)

 {

     int kk,m,l;

     float u[],p,q;

     // 初值

     memset(s, , *sizeof(float));

     // 特例处理

     if (n < )

     {

         return s[];

     }

     if (n == )

     {

         s[] = YfGetFloatValue(valuesPtr, stride, );

         s[] = s[];

         return s[];

     }

     float xStep = 1.0f/(n - );

     if (n == )

     {

         float y0 = YfGetFloatValue(valuesPtr, stride, );

         float y1 = YfGetFloatValue(valuesPtr, stride, );

         s[] = y0;

         s[] = (y1-y0)/xStep;

         s[] = y0 + (y1 - y0)*t;

         return s[];

     }

     // 插值

     if (k<)

     {

         if (t <= xStep)

             kk = ;

         else if (t >= (n-)*xStep)

             kk = n-;

         else

         {

             kk = ;

             m = n;

             while (((kk-m)!=)&&((kk-m)!=-))

             {

                 l=(kk+m)/;

                 if (t < (l-)*xStep)

                     m=l;

                 else

                     kk=l;

             }

             kk=kk-;

         }

     }

     else

         kk=k;

     if (kk>=n-)

         kk=n-;

     u[]=(YfGetFloatValue(valuesPtr, stride, kk+) - YfGetFloatValue(valuesPtr, stride, kk))/xStep;

     if (n==)

     {

         if (kk==)

         {

             u[]=(YfGetFloatValue(valuesPtr, stride, ) - YfGetFloatValue(valuesPtr, stride, ))/xStep;

             u[]=2.0f*u[]-u[];

             u[]=2.0f*u[]-u[];

             u[]=2.0f*u[]-u[];

         }

         else

         {

             u[]=(YfGetFloatValue(valuesPtr, stride, ) - YfGetFloatValue(valuesPtr, stride, ))/xStep;

             u[]=2.0f*u[]-u[];

             u[]=2.0f*u[]-u[];

             u[]=2.0f*u[]-u[];

         }

     }

     else

     {

         if (kk<=)

         {

             u[]=(YfGetFloatValue(valuesPtr, stride, kk+) - YfGetFloatValue(valuesPtr, stride, kk+))/xStep;

             if (kk==)

             {

                 u[]=(YfGetFloatValue(valuesPtr, stride, ) - YfGetFloatValue(valuesPtr, stride, ))/xStep;

                 u[]=2.0f*u[]-u[];

                 if (n==)

                     u[]=2.0f*u[]-u[];

                 else

                     u[]=(YfGetFloatValue(valuesPtr, stride, ) - YfGetFloatValue(valuesPtr, stride, ))/xStep;

             }

             else

             {

                 u[]=2.0f*u[]-u[];

                 u[]=2.0f*u[]-u[];

                 u[]=(YfGetFloatValue(valuesPtr, stride, ) - YfGetFloatValue(valuesPtr, stride, ))/xStep;

             }

         }

         else if (kk>=(n-))

         {

             u[]=(YfGetFloatValue(valuesPtr, stride, kk) - YfGetFloatValue(valuesPtr, stride, kk-))/xStep;

             if (kk==(n-))

             {

                 u[]=(YfGetFloatValue(valuesPtr, stride, n-) - YfGetFloatValue(valuesPtr, stride, n-))/xStep;

                 u[]=2.0f*u[]-u[];

                 if (n==)

                     u[]=2.0f*u[]-u[];

                 else

                     u[]=(YfGetFloatValue(valuesPtr, stride, kk-) - YfGetFloatValue(valuesPtr, stride, kk-))/xStep;

             }

             else

             {

                 u[]=2.0f*u[]-u[];

                 u[]=2.0f*u[]-u[];

                 u[]=(YfGetFloatValue(valuesPtr, stride, kk-) - YfGetFloatValue(valuesPtr, stride, kk-))/xStep;

             }

         }

         else

         {

             u[]=(YfGetFloatValue(valuesPtr, stride, kk  ) - YfGetFloatValue(valuesPtr, stride, kk-))/xStep;

             u[]=(YfGetFloatValue(valuesPtr, stride, kk-) - YfGetFloatValue(valuesPtr, stride, kk-))/xStep;

             u[]=(YfGetFloatValue(valuesPtr, stride, kk+) - YfGetFloatValue(valuesPtr, stride, kk+))/xStep;

             u[]=(YfGetFloatValue(valuesPtr, stride, kk+) - YfGetFloatValue(valuesPtr, stride, kk+))/xStep;

         }

     }

     s[] = fabs(u[]-u[]);

     s[] = fabs(u[]-u[]);

     if ((s[]+1.0f == 1.0f) && (s[]+1.0f == 1.0f))

         p = (u[]+u[])/2.0f;

     else

         p = (s[]*u[]+s[]*u[])/(s[]+s[]);

     s[] = fabs(u[]-u[]);

     s[] = fabs(u[]-u[]);

     if ((s[]+1.0f==1.0f) && (s[]+1.0f==1.0f))

         q = (u[]+u[])/2.0f;

     else

         q = (s[]*u[]+s[]*u[])/(s[]+s[]);

     s[] = YfGetFloatValue(valuesPtr, stride, kk);

     s[] = p;

     s[] = xStep;

     s[] = (3.0f*u[]-2.0f*p-q)/s[];

     s[] = (q+p-2.0f*u[])/(s[]*s[]);

     //if (k<0)

     {

         p=t-(kk*xStep);

         s[]=s[]+s[]*p+s[]*p*p+s[]*p*p*p;

     }

     return s[];

 }

切图：

、

相关软件的下载地址为：http://files.cnblogs.com/WhyEngine/TestSpline.zip