欢迎讨论三次样条曲线的算法！！！

Question

偶最近遇到了一个难题，
关于三次样条曲线的算法问题
输入一系列点，然后绘制出一条三次样条曲线，
欢迎高手来提供源代码予以帮助，
不胜感激！！

Answer 1

//文件名：splineinsert.cpp
/*************************************************************************
*
* 函数名称：
*   splineinsert()
*
* 参数:
* InContourSide        - 存储内轮廓的节点链表
* 返回值:
*   BOOL               - 运算成功返回TRUE，否则返回FALSE。
*
* 说明:
* 该函数用于输出提取出的轮廓数据。
* 
* 要求轮廓图象已经过扫描，建立起轮廓图象，对轮廓数据进行三次样条插值。
************************************************************************///通过此函数调用三次样条插值
BOOL WINAPI splineinsert(ChainHead*  InContourSide,ChainHead* ReplaceSide)
{
int i,j,k;
int direction;
int LongSize;
int TotalSize,other;
ChainHead*  pChain;
Node*       pNode = NULL; 
Node*       mNode = NULL;
Node*       sNode = NULL;
Node*       BreakNode = NULL; int         MidFlag;
Node*       MinNode = NULL;
    Node*       MaxNode = NULL;
//边界条件,即假定边界为切矢无穷大
double  yp1 = 1.7E30;
double  ypn = 1.7E30;

/*
此处代码省,目的是把轮廓分成上下两条
因为样条插值要求沿一个方向是增大/减小的
*/ //x,y 是输入的坐标值 n是数据点的个数 
//yp1,ypn是端始点的二阶导数值 y2为输出的二阶导数值
spline(Xb, Yb, k, yp1, ypn, y2); LongSize = (int)(MaxNode->x-MinNode->x)+1;  //插值点的个数
double*    sig = new double[LongSize];
double*    xcha = new double[LongSize];

xcha[0] = Xb[0];
for(i=1; i<LongSize; i++)
{
xcha[i] = xcha[i-1]+1;
}

for(i=0; i<=LongSize; i++)
{
splint(Xb, Yb, y2, k, xcha[i], sig[i]);
}/*
此处代码省,目的是把轮廓分成上下两条
因为样条插值要求沿一个方向是增大/减小的
*/
//x,y 是输入的坐标值 n是数据点的个数 
//yp1,ypn是端始点的二阶导数值 y3为输出的二阶导数值
spline(Xc, Yc, other, yp1, ypn, y3); for(i=0; i<=LongSize; i++)
{
splint(Xc, Yc, y3, other, xcha[i], sig[i]);
}

for(i=1,j=0; i<LongSize; i++,j++)
{
mNode = new Node[1];
pNode->nextNode = mNode;
pNode = pNode->nextNode;
pNode->x = xcha[j];
pNode->y = sig[j];
}
pNode->nextNode = NULL;
double    ypoint;
Node*  tempNode=NULL; return true;
}//计算出二次导数
void WINAPI spline(double x[], double y[], int n, double yp1,double ypn, double y2[])
{
    double u[100],aaa,sig,p,bbb,ccc,qn,un;
int i,k;
    if( yp1 > 9.9e+29 )
{
        y2[1] = 0;
        u[1] = 0;
}
    else
{
        y2[1] = -0.5;
        aaa = (y[2] - y[1]) / (x[2] - x[1]);
        u[1] = (3.0 / (x[2] - x[1])) * (aaa - yp1);
    }

    for (i = 2; i<=n-1; i++)
{
        sig = (x[i] - x[i - 1]) / (x[i + 1] - x[i - 1]);
        p = sig * y2[i - 1] + 2.0;
        y2[i] = (sig - 1.0) / p;
        aaa = (y[i + 1] - y[i]) / (x[i + 1] - x[i]);
        bbb = (y[i] - y[i - 1]) / (x[i] - x[i - 1]);
        ccc = x[i + 1] - x[i - 1];
        u[i] = (6.0 * (aaa - bbb) / ccc - sig * u[i - 1]) / p;
    }

    if (ypn > 9.9e+29 )
{
        qn = 0.0;
        un = 0.0;
}
    else
{
        qn = 0.5;
        aaa = ypn - (y[n] - y[n - 1]) / (x[n] - x[n - 1]);
        un = (3.0/ (x[n] - x[n - 1])) * aaa;
    }
    y2[n] = (un - qn * u[n - 1]) / (qn * y2[n - 1] + 1.0);

    for (k = n - 1;k>=1;k--)
        y2[k] = y2[k] * y2[k + 1] + u[k];
}
//计算出相应的三次样条插值点
void WINAPI splint(double xa[], double ya[], double y2a[], int n, double& x, double& y)
{
int klo,khi,k;
double h,a,b,aaa,bbb;
    klo = 1;
    khi = n;
loop:   if (khi - klo > 1 )
{
k = (khi + klo) / 2;
if (xa[k] > x)
khi = k;
else
{
klo = k;
}
goto loop;
}
h = xa[khi] - xa[klo];
ASSERT(h!=0); a = (xa[khi] - x) / h;
b = (x - xa[klo]) / h;
aaa = a * ya[klo] + b * ya[khi];
bbb = (a*a*a - a) * y2a[klo] + (b*b*b - b) * y2a[khi];
y = aaa + bbb * h*h /6.0;
}

Answer 2

好复杂阿，有没有一个完整的源程序让我参考一下！！拜托!!!
法到我的邮箱里
[email protected]
谢谢，谢谢！！