C2连续的三次B样条插值（c++）

主要代码如下：
//////////////////////////////////////////////////////////
// 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]);
 }
}
算法实现参考《计算机辅助几何设计》王国瑾 汪国昭 郑建民

本文地址：http://com.8s8s.com/it/it23988.htm