样条之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[5],p,q;

    // 初值
    memset(s, 0, 5*sizeof(float));

    // 特例处理
    if (n < 1) 
    {
        return s[4];
    }
    if (n == 1) 
    { 
        s[4] = YfGetFloatValue(valuesPtr, stride, 0);  
        s[0] = s[4];
        return s[4];
    }

    float xStep = 1.0f/(n - 1);

    if (n == 2)
    { 
        float y0 = YfGetFloatValue(valuesPtr, stride, 0); 
        float y1 = YfGetFloatValue(valuesPtr, stride, 1); 
        s[0] = y0;
        s[1] = (y1-y0)/xStep;
        s[4] = y0 + (y1 - y0)*t;
        return s[4];
    }

    // 插值
    if (k<0)
    { 
        if (t <= xStep) 
            kk = 0;
        else if (t >= (n-1)*xStep) 
            kk = n-2;
        else
        { 
            kk = 1; 
            m = n;
            while (((kk-m)!=1)&&((kk-m)!=-1))
            { 
                l=(kk+m)/2;
                if (t < (l-1)*xStep) 
                    m=l;
                else 
                    kk=l;
            }

            kk=kk-1;
        }
    }
    else 
        kk=k;

    if (kk>=n-1) 
        kk=n-2;

    u[2]=(YfGetFloatValue(valuesPtr, stride, kk+1) - YfGetFloatValue(valuesPtr, stride, kk))/xStep;
    if (n==3)
    { 
        if (kk==0)
        { 
            u[3]=(YfGetFloatValue(valuesPtr, stride, 2) - YfGetFloatValue(valuesPtr, stride, 1))/xStep;
            u[4]=2.0f*u[3]-u[2];
            u[1]=2.0f*u[2]-u[3];
            u[0]=2.0f*u[1]-u[2];
        }
        else
        { 
            u[1]=(YfGetFloatValue(valuesPtr, stride, 1) - YfGetFloatValue(valuesPtr, stride, 0))/xStep;
            u[0]=2.0f*u[1]-u[2];
            u[3]=2.0f*u[2]-u[1];
            u[4]=2.0f*u[3]-u[2];
        }
    }
    else
    { 
        if (kk<=1)
        { 
            u[3]=(YfGetFloatValue(valuesPtr, stride, kk+2) - YfGetFloatValue(valuesPtr, stride, kk+1))/xStep;
            if (kk==1)
            { 
                u[1]=(YfGetFloatValue(valuesPtr, stride, 1) - YfGetFloatValue(valuesPtr, stride, 0))/xStep;
                u[0]=2.0f*u[1]-u[2];
                if (n==4) 
                    u[4]=2.0f*u[3]-u[2];
                else 
                    u[4]=(YfGetFloatValue(valuesPtr, stride, 4) - YfGetFloatValue(valuesPtr, stride, 3))/xStep;
            }
            else
            { 
                u[1]=2.0f*u[2]-u[3];
                u[0]=2.0f*u[1]-u[2];
                u[4]=(YfGetFloatValue(valuesPtr, stride, 3) - YfGetFloatValue(valuesPtr, stride, 2))/xStep;
            }
        }
        else if (kk>=(n-3))
        { 
            u[1]=(YfGetFloatValue(valuesPtr, stride, kk) - YfGetFloatValue(valuesPtr, stride, kk-1))/xStep;
            if (kk==(n-3))
            { 
                u[3]=(YfGetFloatValue(valuesPtr, stride, n-1) - YfGetFloatValue(valuesPtr, stride, n-2))/xStep;
                u[4]=2.0f*u[3]-u[2];
                if (n==4) 
                    u[0]=2.0f*u[1]-u[2];
                else 
                    u[0]=(YfGetFloatValue(valuesPtr, stride, kk-1) - YfGetFloatValue(valuesPtr, stride, kk-2))/xStep;
            }
            else
            { 
                u[3]=2.0f*u[2]-u[1];
                u[4]=2.0f*u[3]-u[2];
                u[0]=(YfGetFloatValue(valuesPtr, stride, kk-1) - YfGetFloatValue(valuesPtr, stride, kk-2))/xStep;
            }
        }
        else
        { 
            u[1]=(YfGetFloatValue(valuesPtr, stride, kk  ) - YfGetFloatValue(valuesPtr, stride, kk-1))/xStep;
            u[0]=(YfGetFloatValue(valuesPtr, stride, kk-1) - YfGetFloatValue(valuesPtr, stride, kk-2))/xStep;
            u[3]=(YfGetFloatValue(valuesPtr, stride, kk+2) - YfGetFloatValue(valuesPtr, stride, kk+1))/xStep;
            u[4]=(YfGetFloatValue(valuesPtr, stride, kk+3) - YfGetFloatValue(valuesPtr, stride, kk+2))/xStep;
        }
    }

    s[0] = fabs(u[3]-u[2]);
    s[1] = fabs(u[0]-u[1]);
    if ((s[0]+1.0f == 1.0f) && (s[1]+1.0f == 1.0f))
        p = (u[1]+u[2])/2.0f;
    else 
        p = (s[0]*u[1]+s[1]*u[2])/(s[0]+s[1]);

    s[0] = fabs(u[3]-u[4]);
    s[1] = fabs(u[2]-u[1]);
    if ((s[0]+1.0f==1.0f) && (s[1]+1.0f==1.0f))
        q = (u[2]+u[3])/2.0f;
    else 
        q = (s[0]*u[2]+s[1]*u[3])/(s[0]+s[1]);

    s[0] = YfGetFloatValue(valuesPtr, stride, kk);
    s[1] = p;
    s[3] = xStep;
    s[2] = (3.0f*u[2]-2.0f*p-q)/s[3];
    s[3] = (q+p-2.0f*u[2])/(s[3]*s[3]);

    //if (k<0)
    { 
        p=t-(kk*xStep);
        s[4]=s[0]+s[1]*p+s[2]*p*p+s[3]*p*p*p;
    }

    return s[4];

}

切图：

 、 

 

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