CAD中的spline详解

从dxf文件中提取点、直线、圆、弧等元素比较简单，但是Spline的处理比较麻烦。经过一段时间探索总结一下成果。

一、基本公式

1.有理样条曲线

查阅一些资料，认为CAD中使用的Spline 是非均匀有理样条曲线。实测CAD中每个控制点权重都是-1，所以下面的有理样条公式蜕变成标准B样条公式。

2.B样条曲线

B样条曲线的表达式如下

d-多项式阶数，CAD中默认阶数4，方程最高次幂degree=3

k-对应控制点数，比如n+1个控制点，则k取值[0,n]

u-参数方程中的“参数”，多个u值构成节点向量，

在CAD中节点向量个数=控制点数+阶数，所以节点向量最后索引是n+d,CAD中的节点向量值是根据拟合点评估的，后面举例说明。

3.混合函数

B样条曲线的混合函数由Cox-deBoor递归公式定义

推导一次导数（没有使用）

推导二次导数（没有使用）

二、CAD样条实例

1.绘制spline

在CAD中使用Spline命令随便绘制一段简单的样条曲线，如下图：

如上图所示，CAD中样条曲线拟合点数量3，阶数4，那么控制点数就是3+4-2。而且默认参数新生成的样条曲线前三个控制点共线，后三个控制点共线。而且线段比例是可以计算出来的，这一点对实现自己的样条非常关键。

2.参数方程

上图是一个3次多项式,阶数d=4，控制点数5，则k取值[0,4]，控制点是C0, C1, C2, C3, C4，将这些参数带入公式，可得该曲线的表达式：

3.计算混合函数

根据上面的公式可以总结出这样一个表：

4.混合函数计算代码


double zmSpline::blend(int k, int d, double u)
{double res = 0;if( d == 1){//1阶时取值范围只有0和1//教科书给的公式   U_k<=u<=U_k+1 ,在节点处会带来问题，//=不能两边同时存在，终点还要特殊处理if(u < m_knots[m_knots.size() - 1]) {res = (m_knots[k] <= u && u < m_knots[k + 1]) ? 1 : 0;}else {res = (m_knots[k] < u && u <= m_knots[k + 1]) ? 1 : 0;}}else{double div1 = m_knots[k + d - 1] - m_knots[k];double div2 = m_knots[k + d] - m_knots[k + 1];//如果分母为0，认为0/0=0double c1 = (std::abs(div1 ) < 10e-14) ? 0 : (u - m_knots[k]) / div1;double c2 = (std::abs(div2) < 10e-14) ? 0 : (m_knots[k + d] - u) / div2;res = c1 * blend(k, d - 1, u) + c2 * blend(k + 1, d - 1, u);}return res;}

5.计算节点向量

如果是单纯的读取dxf文件，节点向量可以从文件读取，是已知量，不用算，下面描述节点向量如何根据拟合点来计算。

控制点数5，阶数4，节点向量个数5+4=9，也就是：

绘制样条时已知量只有这3个拟合点，控制点是不知道的。节点向量也不知道。默认是根据拟合点弦长（还有其他方法，弦长平方根之类的）计算的节点向量。根据上面的弦长可以计算：

上面只得到了3个节点值，但是实际需要9个。前后端点比较特殊，阶数4，那么前后端点对应的节点值需要各自重复4次，所以最终得到的节点向量如下：

刚好是9个，这不是巧合。假设拟合点数量是nFit, 阶数为d,那么控制点数量nCtrl=nFit+(d-2),那么节点数量

nKnot=nCtrl+d= nFit+(d-2)+d= nFit+2*d-2

根据nFit个拟合点可以计算出nFit个节点（第一个是0），前后节点各自重复d次，总的节点数量就是nKnot=nFit+2*d-2。

6.节点向量计算代码


void zmSpline::computeKnotFromFitPoints()
{int size = m_fitPts.size();if(size < 2) {return;}//控制点数要比拟合点数多2//2点拟合-4点控制//5点拟合-7点控制......int n = size + 2;//阶数dint d = m_info.m_degree + 1;//节点数=控制点数+阶数m_knots.resize(n + d, 0);//节点前d个是0for(int i = 0; i < size - 1; i++){double dx = m_fitPts[i + 1].X() - m_fitPts[i].X();double dy = m_fitPts[i + 1].Y() - m_fitPts[i].Y();double length = std::hypot(dx, dy);m_knots[d + i] = m_knots[d + i - 1] + length;}//节点后d个相同，倒数d-1个复制倒数第d个for(int i = m_knots.size() - d + 1; i < m_knots.size(); i++) {m_knots[i] = m_knots[i - 1];}
}

7.读取dxf文件中的Spline

使用libDxf读取dxf文件，能够得到样条曲线的如下信息：

Spline信息结构体代码：

//从DXF文件中获取的信息
typedef struct {unsigned int m_degree;unsigned int m_nKnots;unsigned int m_nControl;unsigned int m_nFit;int m_flags;double m_tangentStartX;double m_tangentStartY;double m_tangentStartZ;double m_tangentEndX;double m_tangentEndY;double m_tangentEndZ;
} SplineInfo;

再看一下这个多项式，所有都是已知量，将u细分带入计算，便能够计算去线上的任意点坐标。u取值范围[0, 135.3347],如果要计算101个曲线点，将u分成100份：0，1.353347，……，分别带入公式计算可以得到101个曲线点，在自己的程序中连接各点便能得到与dxf中一样的样条。当然计算越多，精度越高。

从上面的内容可以实现读取dxf中的spline，然后在自己的代码中显示曲线。如果是要在自己的代码绘制spline，需要继续往下探索。

将u的取值范围划分50份，逐个计算得到50个点，然后将点拷贝到CAD中可以看到计算结果与原始曲线的拟合程度。

三、绘制自己的spline

1.已知量

刚开始绘制的时候一切都是未知的，用鼠标随机生成的点就是spline的拟合点。然后我们默认阶数是4，与CAD中一致。假如我们用鼠标在自己的程序中获取了如下的3个拟合点：

此时我们的已知量如下：

a.阶数d=4

b.拟合点数nFit=3

c.控制点数nCtrl=nFit+d-2=5

d.节点数量nKnots=nCtrl+d=9,还可以用弦长评估出这个9个值

e. 可以计算出来

要求解曲线的参数方程，需要根据上面的已知量，计算出5个控制点。

5个未知数，只有3个方程，还需要2个方程,其实是已知的，可以将5个未知数削减为3个，上面的3个方程只保留中间1个，为了后面描述的完整性，还是按照5个未知数计算。

2.边界条件

查阅一些资料，大多数是通过评估或者设定端点的导数来增加两个方程，如果在绘制spline时，我们手动输入了端点的切向，那么就要通过计算导数来增加两个方程。如果采用默认值，切向是未知的，所以计算导数也不能用。经过一段时间思考和观察，发现CAD中前3个控制点、后3个控制点总是共线（前提是绘制的时候没有输入端点切向）。

观察发现下面的比例关系，于是有了某种猜测：

再随便画一个复杂点的spline观察

再次验证了：前后三个控制点不但共线，而且线段比例和拟合点弦长有关系，而拟合点的弦长又推算出节点向量。所以前后三个控制点的长度比例和可以用节点向量计算。（这是我个人观察和猜测，毕竟没法获取CAD中真正的计算方式，不保证真的正确）

把上面的猜测应用到那条简单的spline里，往下计算：

已知节点向量：

我们得到了下面的方程：

规范一下样式：

再加上前面的3个方程，刚好5个方程，5个未知数，汇总如下：

整理成矩阵样式:

再次说明一下，是已知量，可以削减一下上面的矩阵，但是为了完整描述，仍然把当做未知量计算。

上述过程简化描述如下：

a.默认曲线阶数d

b.获取nFit个拟合点

c.计算nFit+2*d-2个节点值

d.计算nFit+d-2个控制点系数

e.应用前后3个控制点共线、线段比例这个2个边界条件

f.整理成矩阵，计算出nFit+d-2控制点

g.整理出样条参数方程

3.从拟合点计算控制点的代码

void zmSpline::computeCtrlFromFitPoints()
{//dxf中控制点数=拟合点数+2//控制点头尾就是拟合点头尾const int size = m_fitPts.size() + 2;m_ctrlPts.resize(size, SplinePoint());//构建系数矩阵，size个控制点要构建size×size的系数矩阵using namespace Eigen;MatrixXd A(size, size);A.setZero();//所以u的索引是从d-1开始的//i行代表拟合点i对应的计算系数//k列代表控制点k对应的混合函数int d = m_info.m_degree + 1;for(int i = 0; i < size - 2; i++)for(int k = 0; k < size; k++) {A(i, k) = blend(k, d, m_knots[i + d - 1]);}//假如有6个拟合点，就有8个控制点，阶数4，则有8+4=12个节点，//比如0，0，0，0，5，9，15，26，30，30，30，30//控制点C0,C1,C2,C3,C4,C5,C6,C7//初次生成的曲线，前三个控制点共线//(C0-C1)=A0*(C1-C2)   ->   C0+(-1-A0)*C1+A0*C2=0//A0=(5-0)/9//初次生成的曲线，后三个控制点共线//(C7-C6)=A1*(C6-C5)   ->   A1*C5+(-1-A1)*C6+C7=0//A1=(30-26)/(30-15)double A0 = (m_knots[d] - m_knots[d - 1]) / m_knots[d + 1];double A1 = (m_knots[size] - m_knots[size - 1]) / (m_knots[size] - m_knots[size - 2]);A(size - 2, 0) = 1;A(size - 2, 1) = -1 - A0;A(size - 2, 2) = A0;A(size - 1, size - 3) = A1;A(size - 1, size - 2) = -1 - A1;A(size - 1, size - 1) = 1;
//    std::cout << "A:\n" << A << std::endl;//构建拟合点矩阵，最后2行是0MatrixXd B(size, 3);B.setZero();for(int i = 0; i < m_fitPts.size(); i++){B(i, 0) = m_fitPts[i].X();B(i, 1) = m_fitPts[i].Y();B(i, 2) = m_fitPts[i].Z();}/*zmVector v0, v1, v2;BesselTanget(m_fitPts[0], m_fitPts[1], m_fitPts[2], v0, v1, v2);//起点导数B(size - 2, 0) = v0.X() ;B(size - 2, 1) = v0.Y();B(size - 2, 2) = v0.Z();zmVector v_2, v_1, v_0;int sizeOfFit = m_fitPts.size();BesselTanget(m_fitPts[sizeOfFit - 3], m_fitPts[sizeOfFit - 2], m_fitPts[sizeOfFit - 1], v_2, v_1, v_0);//终点导数B(size - 1, 0) = v_0.X();B(size - 1, 1) =  v_0.Y();B(size - 1, 2) = v_0.Z();*///    std::cout << "B:\n" << B << std::endl << std::endl;//求解控制点矩阵MatrixXd X(size, 3);X = A.fullPivLu().solve(B);
//    std::cout << "X:\n" << X << std::endl << std::endl;for(int i = 0; i < size; i++){double x = X(i, 0);double y = X(i, 1);double z = X(i, 2);SplinePoint point(x, y, z, -1);m_ctrlPts[i] = point;}
}

计算出控制点之后带入B样条曲线，便能计算曲线上各个点。