C++实现Rhino中画准均匀B样条曲线功能

1. B样条曲线特征

B样条曲线是由Bezier曲线演变而来的，其最大的优点是通过调整局部控制点来改变曲线形状。对于n次B样条曲线，其公式如下。可以看出，对于3次即4阶B样条曲线，至少需要4个控制点。
对于下面的公式，我们已知的是: 所有控制点 P_i 和曲线阶数n。需要求解的有: 每个控制点对应的B_i,n和节点向量u。

2. 非均匀B样条有理曲线

对于B样条曲线的基础知识，可以参考其他文章，本文章是默认对B样条曲线有一定的了解，目的是用C++实现Rhino中非均匀B样条曲线的绘制。
对于非均匀B样条曲线，其节点向量间隔为定值，为了方便计算，此处设置节点向量间隔为1。对于k阶非均匀B样条曲线，节点向量起始为0，重复k次，末端节点也重复k次。节点向量计算代码如下：

例如：对于3次4阶曲线，控制点为5个，控制点坐标如下:

序号	坐标值
1	(24.3785, -2.87838, 0)
2	(4.82241, -9.35762, 12.9518)
3	(22.7527, 17.7923, 0)
4	(-4.69538, 16.966, -9.8671)
5	(-10.6367, -10.1958, 0)

其节点向量计算结果为：
{0,0,0,0,1,2,2,2,2}

在Rhino中，其曲线如下图所示：

//k为曲线阶数，point_num为控制点数量，返回为节点向量
vector<int> calNodes(const int& k, const int& point_num)
{
  vector<int> nodes;
  for (int i = 0; i < k - 1; i++)
    nodes.push_back(0);
  for (int i = 0; i < k + point_num - 2 * k + 2; i++)
    nodes.push_back(i);
  for (int i = 0; i < k - 1; i++)
    nodes.push_back(nodes.back());
  return nodes;
}

3. 计算每个控制点对应的系数

对于曲线上的每一点，都会对应一个u值，寻找其u在节点向量中的位置，根据**Cox-deBoor**递推公式，迭代计算其对应的系数B_i,k，代码如下：

其对应的代码如下，根据控制点序号，曲线阶数，曲线上点的位置和节点向量，可以计算处这个控制点对应的系数矩阵。

  //i为对应的控制点序号，k为曲线阶数，u为曲线上的位置点，nodes为曲线节点向量
Vector4d  calCoffe(const int& i, const int& k, const double& u, const vector<int>& nodes)
{
  Vector4d coffe;
  if (1 == k)
  {
    if ((u >= nodes[i] && u < nodes[i + 1]) || fabs(nodes.back() - u) < 1e-8)
      coffe = Vector4d(0, 0, 0, 1);
    else
      coffe = Vector4d::Zero();
  }
  else
  {
    double length1 = nodes[i + k - 1] - nodes[i];
    double length2 = nodes[i + k] - nodes[i + 1];
    if (0.0 == length1)
      length1 = 1.0;
    if (0.0 == length2)
      length2 = 1.0;
    Vector4d tmp1, tmp2;
    tmp1(0) = -calCoffe(i, k - 1, u, nodes)(0) * nodes[i] / length1 + calCoffe(i, k - 1, u, nodes)(1) / length1;
    tmp1(1) = -calCoffe(i, k - 1, u, nodes)(1) * nodes[i] / length1 + calCoffe(i, k - 1, u, nodes)(2) / length1;
    tmp1(2) = -calCoffe(i, k - 1, u, nodes)(2) * nodes[i] / length1 + calCoffe(i, k - 1, u, nodes)(3) / length1;
    tmp1(3) = -calCoffe(i, k - 1, u, nodes)(3) * nodes[i] / length1;

    tmp2(0) = calCoffe(i + 1, k - 1, u, nodes)(0) * nodes[i + k] / length2 - calCoffe(i + 1, k - 1, u, nodes)(1) / length2;
    tmp2(1) = calCoffe(i + 1, k - 1, u, nodes)(1) * nodes[i + k] / length2 - calCoffe(i + 1, k - 1, u, nodes)(2) / length2;
    tmp2(2) = calCoffe(i + 1, k - 1, u, nodes)(2) * nodes[i + k] / length2 - calCoffe(i + 1, k - 1, u, nodes)(3) / length2;
    tmp2(3) = calCoffe(i + 1, k - 1, u, nodes)(3) * nodes[i + k] / length2;
    coffe = tmp1 + tmp2;
  }
  return coffe;
}

4. 计算曲线上的点

最后，我们来创建这条3次4阶，拥有5个控制点的非均匀有理B样条曲线，在曲线上均匀采样200个点。其计算步骤如下：

1. 首先通过calNodes函数计算其节点向量；
2. 针对曲线上节点坐标为u的点，通过calCoffe函数计算对应的控制点系数矩阵；
3. 将每个控制点和对应的系数矩阵相乘，得到曲线上节点坐标为u时的空间坐标point；
4. 重复2和3的操作，直至采样阶数，共计200个点.

#include <iostream>
#include <vector>
#include <Eigen\Dense>
using namespace std;
using namespace Eigen;
int main(int argc, const char* argv[])
{
  vector<Vector3d> cv_points;                               //曲线控制点
  cv_points.push_back(Vector3d(24.3785, -2.87838, 0));      //第1个控制点
  cv_points.push_back(Vector3d(4.82241, -9.35762, 12.9518));//第2个控制点
  cv_points.push_back(Vector3d(22.7527, 17.7923, 0));       //第3个控制点
  cv_points.push_back(Vector3d(-4.69538, 16.966, -9.8671)); //第4个控制点
  cv_points.push_back(Vector3d(-10.6367, -10.1958, 0));     //第5个控制点

  int k = 4;                                                //设置曲线阶数为4，曲线最高次数为3，共5个控制点
  vector<int> nodes = calNodes(k, cv_points.size());        //计算曲线节点向量结果为{0,0,0,1,2,3,3,3}

  int range[2] = { nodes[k - 1],nodes[cv_points.size()] };  //得到曲线上u的取值范围为0-3
  vector<Vector3d> points;                                  //用来存储曲线上的坐标点
  for (double u = range[0]; u < range[1]; u += 0.01)       //求曲线上的点，采样间隔为0.01，共得到300个点
  {
    Vector4d U(u * u * u, u * u, u, 1);
    Vector3d point(0, 0, 0);                                //存储节点为u时的坐标
    int num = floor(u - range[0]);
    if (u == range[1])
      num--;
    for (int i = num; i < num + k; i++)
    {
      Vector4d coffe = calCoffe(i, k, u, nodes);
      point += coffe.dot(U) * cv_points[i];
    }
    points.push_back(point);
    cout << point(0) << "," << point(1) << "," << point(2) << endl;
  }
}

将200个点画入Rhino中，点落在曲线上，所以计算正确。

5.总结

本文主要针对简单的非均匀有理B样条曲线，如有相关深入的问题，欢迎探讨。后续还会有进阶版的关于B样条曲线的更新。