样条之连分式插值函数

核心代码：

//////////////////////////////////////////////////////////////////////
// 连分式等距插值
//////////////////////////////////////////////////////////////////////
static float GetValuePqs(const void* valuesPtr, int stride, int n, float t)
{ 
    int i,j,k,m, l;
    float z,hh,xi,xj;
    float b[8];

    // 初值
    z = 0.0;

    // 特例处理
    if (n < 1) 
    {
        return(z);
    }
    if (n == 1) 
    { 
        z = YfGetFloatValue(valuesPtr, stride, 0);  
        return(z);
    }
    if (n == 2)
    { 
        float y0 = YfGetFloatValue(valuesPtr, stride, 0); 
        float y1 = YfGetFloatValue(valuesPtr, stride, 1); 
        z = y0 + (y1 - y0)*t;
        return(z);
    }

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

    // 连分式插值
    if (n <= 8) 
    { 
        k = 0; 
        m = n;
    }
    else if (t < (4.0f*xStep)) 
    { 
        k=0; 
        m=8;
    }
    else if (t > ((n-5)*xStep)) 
    { 
        k = n-8; 
        m = 8;
    }
    else 
    { 
        k = (int)(t/xStep)-3; 
        m = 8;
    }

    b[0] = YfGetFloatValue(valuesPtr, stride, k); 
    for (i = 2; i <= m; i++)
    { 
        hh = YfGetFloatValue(valuesPtr, stride, i+k-1); 
        l = 0; 
        j = 1;

        while ((l == 0) && (j <= i-1))
        { 
            if (fabs(hh-b[j-1])+1.0f == 1.0f) 
            {
                l = 1;
            }
            else
            { 
                xi = (i+k-1)*xStep;
                xj = (j+k-1)*xStep;
                hh = (xi-xj)/(hh-b[j-1]);
            }

            j = j+1;
        }

        b[i-1]=hh;
        if (l != 0) 
        {
            b[i-1] = 1.0e+35F;
        }
    }

    z = b[m-1];
    for (i = m-1; i >= 1; i--)
    {
        z = b[i-1]+(t-(i+k-1)*xStep)/z;
    }

    return(z);
}

切图：

 

 

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