主要代码如下：
//////////////////////////////////////////////////////////
// cubic B-spline Interpolate
// n : = num-1,means n segemnt of the knot-zone
// made at 01/11/2004 by jzp
//////////////////////////////////////////////////////////
#define INTERPDGR 3
#define _DelPointGroup_(x) { assert((x)); delete [](x); (x)=NULL; }
void CSpline::InterpolateValue(int n)
{
CPoint3D * InterPoints;
InterPoints = new CPoint3D [n+3];
// convey
Convey2InterpolatePoints();
//
InterPoints[0] = pInterpolatePoints[0];
for (int i = 2; i <= n; i ++)
{
InterPoints[i] = pInterpolatePoints[i-1];
}
InterPoints[n+2] = pInterpolatePoints[n];
// assign the B-Spline object
SetDegree(INTERPDGR);
mKnotNumber = n+1;
SetCPNumber(n+INTERPDGR);
AssignMemory();
// define the coefficient matrix
double f3,g3,coef,temp,*TriDiaMatrix[3];
for (i = 0; i < 3; i ++)
{
TriDiaMatrix[i] = new double [n+3];
}
//_________________________________
// assign the knot value
double edgelength;
CPoint3D Edge;
double *t;
t = new double[n+1];
t[0] = 0.0;
edgelength = SplineEdgeLong(n+1);
/*for (i = 1; i <= n; i ++)
{
Edge = pInterpolatePoints[i] - pInterpolatePoints[i-1];
edgelength += Edge.length();
}*/
for (i = 1; i <= n; i ++)
{
Edge = pInterpolatePoints[i] - pInterpolatePoints[i-1];
t[i] = t[i-1] + Edge.length()/edgelength;
}
// assign the knot vector
pKnotValue[0] = 0.0;
pMultiDegree[0] = INTERPDGR+1;
for (i = 1; i < n; i ++)
{
pKnotValue[i] = t[i];
pMultiDegree[i] = 1;
}
pKnotValue[n] = t[n];
pMultiDegree[n] = INTERPDGR+1;
// assign the coefficient matrix
// \ \ // 0 1 2
// \ \ // swap the two rows at head and tail
TriDiaMatrix[1][0] = 1.0;//CalSplineBasis(0,3,t[0]);
TriDiaMatrix[2][0] = 0.0;
TriDiaMatrix[2][1] = t[1]-t[0];
TriDiaMatrix[0][0] = t[2]-t[0];
TriDiaMatrix[1][1] = 2*t[0]-t[1]-t[2];
TriDiaMatrix[1][n+2] = 1.0;//CalSplineBasis(n+2,3,t[n]);
TriDiaMatrix[0][n+1] = 0.0;
TriDiaMatrix[0][n] = t[n]-t[n-1];
TriDiaMatrix[2][n+1] = t[n]-t[n-2];
TriDiaMatrix[1][n+1] = 2*t[n-2]-t[n-1]-t[n];
// assign the rest rows
for (i = 1; i <= n-1; i ++)
{
TriDiaMatrix[0][i] = CalSplineBasis(i,3,t[i]);
TriDiaMatrix[1][i+1] = CalSplineBasis(i+1,3,t[i]);
TriDiaMatrix[2][i+1] = CalSplineBasis(i+2,3,t[i]);
}
// simplize the coefficient matrix,diagonalize the matrix
for (i = 0; i <= n ; i ++)
{
coef = TriDiaMatrix[0][i]/TriDiaMatrix[1][i];
TriDiaMatrix[0][i] = 0.0;
TriDiaMatrix[1][i+1] -= coef*TriDiaMatrix[2][i];
InterPoints[i+1] = InterPoints[i+1] - InterPoints[i]*coef ;
}
//__________________________________________________
CPoint3D *CtrlPoints;
CtrlPoints = new CPoint3D [n+3];
// calculate the Control points
CtrlPoints[n+2] = InterPoints[n+2]*(1.0/TriDiaMatrix[1][n+2]);
for (i = n+1; i >= 0; i --)
{
CtrlPoints[i] = (InterPoints[i] - CtrlPoints[i+1]*TriDiaMatrix[2][i])*(1.0/TriDiaMatrix[1][i]);
}
memcpy(pControlPoints,CtrlPoints,sizeof(CPoint3D)*(n+3));
_DelPointGroup_(t);
_DelPointGroup_(CtrlPoints);
_DelPointGroup_(InterPoints);
for (i = 0; i < 3; i ++)
{
_DelPointGroup_(TriDiaMatrix[i]);
}
}
算法实现参考《计算机辅助几何设计》王国瑾 汪国昭 郑建民
