CAD中的spline详解

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=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样条曲线，便能计算曲线上各个点。