从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=0
double 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;
//阶数d
int d = m_info.m_degree + 1;
//节点数=控制点数+阶数
m_knots.resize(n + d, 0);
//节点前d个是0
for(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行是0
MatrixXd 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样条曲线，便能计算曲线上各个点。
扫一扫
4707
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
