鞋业CAD软件柔性内核技术:源码[13-

void BSpline::Inverse (int m, Vector *pv)
//功能: 用"追赶法"解三对角方程组,求出控制点.
//      由于第一个和最后一个控制点分别与第一个和最后一个型值点重合,
//      所以必须省略两个方程.
//参数: m -- 点数; pv -- 点
//返回: 无
{
    int     i;
    float   *a, *b, *c, *u, l;
    Vector  *d, *p;

    //如果只有两个型值点,则认为是一线段
    if (m == 2) {
        ctl[0] = pv[0];
        ctl[1] = (pv[1] - pv[0])     / 3 + pv[0];
        ctl[2] = (pv[1] - pv[0]) * 2 / 3 + pv[0];
        ctl[3] = pv[1];
        return;
        }

    a = new float[m + 1];
    b = new float[m + 1];
    c = new float[m + 1];
    u = new float[m + 1];
    d = new Vector[m + 1];
    p = new Vector[m + 1];

    //由B样条方程可推出: p[0]' = 3 (v[1] - v[0]);
    //而 v[0] = p[0];
    //所以 3v[1] = 3v[0] + p[0]'
    //同理 3v[num] = 3v[num+1] - p[num-1]'
    d[1] = EndCondition (pv[0], pv[1], pv[2]) + pv[0] * 3;
    d[m] = EndCondition (pv[m - 1], pv[m - 2], pv[m - 3]) + pv[m - 1] * 3;

    for (i = 2; i < m; i++) d[i] = pv[i - 1];

/*
    |  3    0                                   | | ctl[1]   |   | pv[0]'+3pv[0]       |
    | a(2) b(2) c(2)                            | | ctl[2]   |   | pv[1]               |
    |      a(3) b(3) c(3)                       | | ctl[3]   |   | pv[2]               |
    |           ......                          | |  ...     | = |  ...                |
    |             a(m-2) b(m-2) c(m-2)          | | ctl[m-2] |   | pv[m-3]             |
    |                    a(m-1) b(m-1) c(m-1)   | | ctl[m-1] |   | pv[m-2]             |
    |                             0      3      | | ctl[m]   |   | pv[m-1]'+3pv[num-1] |
*/
    a[1] = 0.0; b[1] = 3.0; c[1] = 0.0;
    for (i = 1; i < m - 1; i++) {
        a[i + 1] = coe[i][0][0];
        b[i + 1] = coe[i][1][0];
        c[i + 1] = coe[i][2][0];
        }
    a[m] = 0.0; b[m] = 3.0; c[m] = 0.0;

 //try
 //{
  u[1] = b[1];
  p[1] = d[1];
  for (i = 2; i <= m; i++)
  {
   VERIFY (u[i - 1] != 0);
   l = a[i] / u[i - 1];
   u[i] = b[i] - l * c[i - 1];
   p[i] = d[i] - p[i - 1] * l;
  }

  ctl[0] = pv[0];
  ctl[m + 1] = pv[m - 1];
  ctl[m] = p[m] / u[m];
  for (i = m - 1; i > 0; i--)
  {
   VERIFY (u[i] != 0);
   ctl[i] = (p[i] - ctl[i + 1] * c[i]) / u[i];
  }
 //} catch(...) {}

    delete a;
    delete b;
    delete c;
    delete u;
    delete d;
    delete p;
}

int BSpline::Make (int m, Vector *p, int autodelete)
//功能: 从num个型值点插值出一条B样条曲线
//参数: m -- 点数; p -- 点
//返回: 若有效点数少于3返回false,否则返回true
{
    Flush ();
    if (m < 2 || (m==2&&p[0]==p[1])) {
        if(autodelete) delete p;
        return FALSE;
       }

    t = new float[m + 7];       // increase space for extend
    ōrg = new Vector[m];
    ctl = new Vector[m + 2];
    buf = new Vector[2];        //temporary allocate
    coe = new Coeff[m - 1];
    pos = new int[m];
 posbackup = new int[m];

    num = Knots (m, p, knotlen);
    if (num < 2) {
        Flush ();
        return FALSE;
        }

    //计算系数矩阵
    for (int i = 0; i < num - 1; i++)
  CalCoeff (coe[i], i);

    Inverse (num, org);

    segnum = Split (fit_error_disp, &buf);
 memcpy(posbackup,pos, m*sizeof(int));

    if (autodelete) delete p;
 
 GetMinMax();

    return TRUE;
}
 
Vector BSpline::P (int i, float u) const
//功能: 计算B样条曲线中某段上的一点
//参数: i -- 曲线段, u -- 给定的u值
//返回: 该点坐标
{
    ASSERT (i >= 0 && i < num - 1);
    Vector v =
        ctl[i    ] * (coe[i][0][0] * (1 - u) * (1 - u) * (1 - u)) +
        ctl[i + 1] * (coe[i][1][0] + u * (coe[i][1][1] + u * (coe[i][1][2] + u * coe[i][1][3])))

+
        ctl[i + 2] * (coe[i][2][0] + u * (coe[i][2][1] + u * (coe[i][2][2] + u * coe[i][2][3])))

+
        ctl[i + 3] * (coe[i][3][3] * u * u * u);

    return v;
}

Vector BSpline::TP (int i, float u) const
//功能: 计算B样条曲线中某段上的一点的切矢(导数)
//参数: i -- 曲线段, u -- 给定的u值
//返回: 切矢 (Tangent Vector)
{
    ASSERT (i >= 0 && i < num - 1);
    Vector v =
        ctl[i    ] * (coe[i][0][0] * (1 - u) * (1 - u) * (-3)) +
        ctl[i + 1] * (coe[i][1][1] + u * (coe[i][1][2] * 2 + u * coe[i][1][3] * 3)) +
        ctl[i + 2] * (coe[i][2][1] + u * (coe[i][2][2] * 2 + u * coe[i][2][3] * 3)) +
        ctl[i + 3] * (coe[i][3][3] * u * u * 3);
    //注意: 只有v/(t[i+4]-t[i+3])后才能保证TP(i,1)=TP(i+1,0)
 VERIFY ( (t[i + 4] - t[i + 3]) != 0);
    v /= (t[i + 4] - t[i + 3]);

    return v;
}

Vector BSpline::T2P (int i, float u) const
//功能: 计算B样条曲线中某段上的一点的二阶导数
//参数: i -- 曲线段, u -- 给定的u值
//返回: 二阶导数
{
    ASSERT (i >= 0 && i < num - 1);
    Vector v =
        ctl[i    ] * (coe[i][0][0] * (1 - u) * 6) +
        ctl[i + 1] * (coe[i][1][2] * 2 + u * coe[i][1][3] * 6) +
        ctl[i + 2] * (coe[i][2][2] * 2 + u * coe[i][2][3] * 6) +
        ctl[i + 3] * (coe[i][3][3] * u * 6);
    //注意: 只有v/((t[i+4]-t[i+3])(t[i+4]-t[i+3]))后才能保证T2P(i,1)=T2P(i+1,0)
 VERIFY ( (t[i + 4] - t[i + 3]) * (t[i + 4] - t[i + 3]) != 0);
    v /= (t[i + 4] - t[i + 3]) * (t[i + 4] - t[i + 3]);

    return v;
}

Vector BSpline::NP (int i, float u) const
//功能: 计算B样条曲线中某段上的一点的法矢
//参数: i -- 曲线段, u -- 给定的u值
//返回: 法矢 (Normal Vector)
{
    Vector v = TP (i, u);
    return Vector (v.y, -v.x);
}

void BSpline::ToBezier (int i, Vector p[4]) const
//功能: 将一段B样条曲线转换为Bezier曲线, 解下列线性方程组
//     |  1  0  0  0 |   |P0|   | c00 c10 c20 c30 |   |C0|
//     | -3  3  0  0 | * |P1| = | c01 c11 c21 c31 | * |C1|
//     |  3 -6  3  0 |   |P2|   | c02 c12 c22 c32 |   |C2|
//     | -1  3 -3  1 |   |P3|   | c03 c13 c23 c33 |   |C3|
//     其中, P0, P1, P2, P3是Bezier曲线的控制顶点;
//          C0, C1, C2, C3是B样条曲线的控制顶点
//参数: i -- 曲线段, p -- Bezier曲线的控制顶点
//返回: (无)
{
    int     m, n;
    Vector  d[4];
   
    for (m = 0; m < 4; m++)
        for (n = 0; n < 4; n++)
            d[m] += ctl[i + n] * coe[i][n][m];

    p[0] = d[0];
    p[1] = d[1] / 3 + d[0];
    p[2] = d[2] / 3 + d[1] * 2 / 3 + d[0];
    p[3] = d[3] + d[2] + d[1] + d[0];
}