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

#include "stdafx.h"
#include "bspline.h"
#include "global.h"

#ifdef _DEBUG
#define new DEBUG_NEW
#undef THIS_FILE
static char THIS_FILE[] = __FILE__;
#endif

/*-- Calculate Core -------------------------------------*/

int BSpline::Split (float error, Vector **p) const
//功能: 将每段曲线分为若干段，用于绘图和求交
//参数: error -- 逼近误差
//返回: 点数
{
    #define BUF_LEN 2048

    int     i, k;
    float   a, b;
    Vector  *tmp;
    Vector  n0, n1, p0, p1;
 VectorRun dp;

 if(error <= 0.0025) error = 0.0025f;
 
    //ASSERT (error > 0);
   
    if (previewmode)
 {
        tmp = new Vector[num];
  if(NULL == tmp)
   return -1;

        if (p[0])
   delete p[0];
  p[0] = tmp;

        for (i = 0; i < num; i++)
            p[0][i] = org[i];
        return num;
 }

    tmp = new Vector[BUF_LEN]; //为分割的曲线分配内存
 if(NULL == tmp)
  return -1;

    tmp[0] = P (0, 0); //分割后，曲线段的第一个点
    for (k = 1, i = 0; i < num - 1 && k {
        a = 0;
        b = 1;
        p1 = P (i, 0);  //计算B样条曲线中某段上的一点
        n1 = NP (i, 0); //计算B样条曲线中某段上的一点的法矢

        //计算曲线上某一段的坐标,每段曲线分解为许多段，以ERROR为依据,渐次逼近
        while (a < b && k < BUF_LEN)
  {
            p0 = p1;
            n0 = n1;
            p1 = P (i, b);
            n1 = NP (i, b);

            //采用逼近法计算误差控制范围内的P1点坐标，当计算出的误差小于ERROR时，给出下一个点
            //计算对于曲线段i从p0点到p1点用直线段逼近的误差值,
               while (dp.Delta (p0, n0, p1, n1) > error)
   {
                b = (a + b) / 2;
                p1 = P (i, b);
                n1 = NP (i, b);
            }

            tmp[k++] = p1;
            a = b;
            b = 1;
        }
        pos[i] = k - 1;
    }
    pos[num - 1] = k;

 Vector* tmp2 = NULL;
    tmp2 = new Vector[k];
 if(NULL == tmp2)
 {
  delete [] tmp;
  return -1;
 }

    if (p[0])
  delete p[0];
 p[0] = tmp2;

    memcpy (p[0], tmp, k * sizeof (Vector));

    delete [] tmp;

    return k;
}

int BSpline::SplitNum (float error) const
//功能: 将每段曲线分为若干段，求点数
//参数: error -- 逼近误差
//返回: 点数
{
//    #define BUF_LEN 2048 //2000

    int     i, k;
    float   a, b;
    Vector  n0, n1, p0, p1;
 VectorRun dp;

 if(error <= 0.001) error = 0.001f;
   
    for (k = 1, i = 0; i < num - 1 && k {
        a = 0;
        b = 1;
        p1 = P (i, 0);
        n1 = NP (i, 0);
        while (a < b && k  {
            p0 = p1;
            n0 = n1;
            p1 = P (i, b);
            n1 = NP (i, b);

            while (dp.Delta (p0, n0, p1, n1) > error)
   {
                b = (a + b) / 2;
                p1 = P (i, b);
                n1 = NP (i, b);
            }

            a = b;
            b = 1;
        }
    }

    return k;
}

int BSpline::Knots (int m, Vector *p, float knotlen)
//功能: 计算节点(knots)值t[i]
//      t[i] = 0, 0, 0, 0, l(1), l(2), ... , l(num-1), l(num), l(num), l(num), l(num)
//参数: m -- 点数; p -- 点; knotlen 节点间距,控制曲线光顺度
//返回: 有效型值点数
{
    int     i, j, k;
    float   l=0;
 VectorRun dp;

    t[0] = t[1] = t[2] = t[3] = 0;
    org[0] = p[0];
    for (i = 1, j = 4, k = 1; i < m; i++) {
#ifdef  UNIFORM_BSPLINE
        l += 1;                          //uniform B-spline
#else
        l += dp.Length (p[i] - p[i - 1]);   //non-uniform B-spline ????
#endif
        if (l >= knotlen ) { 
            t[j] = t[j - 1] + l;
            org[k++] = p[i];
            j++;
   l = 0;
  } //99-06-14
  else if(i==m-1) //最后一点
  {
   if(m>2 && k>1) {
    t[j-1] += l;
    org[k-1] = p[i];
   }
   else if(l>=0.0)
   {
    t[j] = t[j - 1] + l;
    org[k++] = p[i];
    j++;
    l = 0;
   }
  }
  //否则l不清零
 }
    t[j] = t[j + 1] = t[j + 2] = t[j - 1];

 if(k==2 && t[4]<=0.0)
  return 0;
    return k;
}

float BSpline::B (int j, int k, float x) const
//功能: 计算x点处B样条基函数的值(利用"de Boor-Cox"递推公式)
//参数: j -- 参数区间([0, num+1])
//      k -- k-1次B样条(此处为4)
//      x -- 变量
//返回: 函数值
//注意: B(j,k,x)的值与coeff是对应的.
//      其中x=t[j+3]处的函数值与coeff[j]相对应:
//          B(j,  4,t[j+3]) == coeff[j][0][0]
//          B(j+1,4,t[j+3]) == coeff[j][1][0]
//          B(j+2,4,t[j+3]) == coeff[j][2][0]
{
    if (k == 1)
        return (x >= t[j] && x < t[j + 1]) ? 1 : 0;
    else {
        float a, b, c, d;

        a = B (j, k - 1, x) * (x - t[j]);
        b = t[j + k - 1] - t[j];
        c = B (j + 1, k - 1, x) * (t[j + k] - x);
        d = t[j + k] - t[j + 1];

  //try
  //{
   //if (a == 0)
   // return (c == 0) ? 0 : c / d;
   //else
   // return (c == 0) ? a / b : a / b + c / d;
   VERIFY (b!=0 && d!=0);
   return (c == 0) ? a / b : a / b + c / d;
  //}
  //catch(...)
  //{
  //}
 }
}

float BSpline::dB (int j, int k, float x) const
//功能: 计算x点处B样条基函数的导数值(利用"de Boor-Cox"递推公式)
//参数: (同上)
//返回: 函数值
{
    float a, b, c, d;

    a = B (j, k - 1, x) * (k - 1);
    b = t[j + k - 1] - t[j];
    c = B (j + 1, k - 1, x) * (k - 1);
    d = t[j + k] - t[j + 1];

 //try
 //{
  //if (a == 0)
  // return (c == 0) ? 0 : c / d;
  //else
  // return (c == 0) ? a / b : a / b + c / d;
 VERIFY (b!=0 && d!=0);
 return (c == 0) ? a / b : a / b + c / d;
 //}
 //catch(...)
 //{
 //}
}

void BSpline::CalCoeff (Coeff& c, int i)
//功能: 计算系数矩阵(4x4)
//      | c[0][0] c[1][0] c[2][0] c[3][0] |
//      | c[0][1] c[1][1] c[2][1] c[3][1] |
//      | c[0][2] c[1][2] c[2][2] c[3][2] |
//      | c[0][3] c[1][3] c[2][3] c[3][3] |
//参数: c -- 系数矩阵, i -- 参数区间([0, num+1])
//返回: 无
{
    float t32, t41, t42, t43, t52, t53, t63;

    t32 = t[i + 3] - t[i + 2];

    t41 = t[i + 4] - t[i + 1];
    t42 = t[i + 4] - t[i + 2];
    t43 = t[i + 4] - t[i + 3];

    t52 = t[i + 5] - t[i + 2];
    t53 = t[i + 5] - t[i + 3];

    t63 = t[i + 6] - t[i + 3];

 //try
 //{
  VERIFY (t41 != 0);
  VERIFY (t42 != 0);
  VERIFY (t52 != 0);
  VERIFY (t53 != 0);
  VERIFY (t63 != 0);

  c[0][0] = t43 * t43 / (t42 * t41);
  c[0][1] = -3 * c[0][0];
  c[0][2] =  3 * c[0][0];
  c[0][3] = -c[0][0];

  c[2][0] = t32 * t32 / (t52 * t42);
  c[2][1] = 3 * t32 * t43 / (t52 * t42);
  c[2][2] = 3 * t43 * t43 / (t52 * t42);
  c[2][3] = - t43 * t43 * (1 / t42 / t52 + 1 / t53 / t63 + 1 / t53 / t52);

  c[3][0] = c[3][1] = c[3][2] = 0;
  c[3][3] = t43 * t43 / (t63 * t53);

  c[1][0] =  1 - c[0][0] - c[2][0];
  c[1][1] =  3 * c[0][0] - c[2][1];
  c[1][2] = -3 * c[0][0] - c[2][2];
  c[1][3] = c[0][0] - c[2][3] - c[3][3];
 //}
 //catch(...)
 //{
 //}
}

Vector BSpline::EndCondition (Vector& p1, Vector& p2, Vector& p3) const
//功能: 计算始末端的端点条件.
//      Bessel(抛物线)方法:
//          认为p1, p2, p3三点构成一抛物线(各点的二阶导数为常数)
//          p1' + p2' = 2 (p2 - p1) / (t2 - t1)
//          p1' + p3' = 2 (p3 - p1) / (t3 - t1)
//          p2' + p3' = 2 (p3 - p2) / (t3 - t2)
//      ==> p1' = (p2 - p1) / (t2 - t1) + (p3 - p1) / (t3 - t1) - (p3 - p2) / (t3 - t2)
//参数: p1 -- 端点;
//      p2 -- 端点后(前)第一点;
//      p3 -- 端点后(前)第二点;
//返回: p  -- 切矢(对于末端是反向的切矢)
{
    float   t21, t32, t31;
 VectorRun dp;

    t21 = dp.Length (p2 - p1);
    t32 = dp.Length (p3 - p2);
    t31 = t21 + t32;

 Vector p;
 //try
 //{
  VERIFY (t21 != 0);
  VERIFY (t31 != 0);
  VERIFY (t32 != 0);
  p = (p2 - p1) / t21 + (p3 - p1) / t31 - (p3 - p2) / t32;
 //}
 //catch(...) {}
    return dp.UnitVector (p) * t21;
}