Simpson积分方法计算NURBS曲线弧长，详细原理+代码实现

本文链接：https://blog.youkuaiyun.com/weixin_45650767/article/details/130140848

文章介绍了使用Simpson积分方法计算NURBS曲线弧长的原理，包括将曲线分段并应用Simpson公式进行数值积分，以及如何处理NURBS曲线的导数。还提供了一个自适应Simpson积分的C++代码示例，用于递归计算满足精度要求的弧长。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

Simpson积分方法计算NURBS曲线弧长，详细原理+代码实现

Simpson 积分方法是一种数值积分方法，可以用于计算曲线的弧长。它的基本思想是将曲线分成若干小段，对每一小段采用 Simpson 公式进行数值积分，然后将各小段的积分结果加起来得到整个曲线的积分结果。
具体地，设曲线

C

由参数方程

\mathbf{r}(t)=(x(t), y(t))

给出，

\in [a, b]

，

n

是分段数，则将

t

分成

n

段，即

t_0=a, t_1, \ldots, t_n=b

，每段长度为

h = (b - a) / n

，则 Simpson 公式可以表示为：

\int_{t_i}^{t_{i+2}} |\mathbf{r}'(t)| \mathrm{d}t \approx \frac{h}{3} (|\mathbf{r}'(t_i)| + 4|\mathbf{r}'(t_{i+1})| + |\mathbf{r}'(t_{i+2})|) \quad (i=0, 2, 4, \ldots, n-2)

这个公式是对

t_i, t_{i+2}]

区间上的函数

|\mathbf{r}'(t)|

进行数值积分的，近似结果可以看作是这个区间上函数曲线下面的面积，即弧长的近似值。将所有区间的积分结果加起来即得到整个曲线的积分结果。
对于 NURBS 曲线，我们可以计算其导数

\mathbf{r}'(t)

，然后代入 Simpson 公式进行数值积分。在计算导数时，可以使用 evaluate_single 函数来计算曲线上某一点的坐标和一阶导数。
下面给出一个示例代码，用自适应 Simpson 积分方法计算 NURBS 曲线的弧长。自适应 Simpson 积分方法是一种递归方法，将整个曲线分成若干小段，对每一小段采用 Simpson 公式进行数值积分，若相邻两次积分的结果相差较大，则将这一小段再分成两半，递归计算，直到精度满足要求为止。代码如下：

#include <iostream>
#include <cmath>
#include <vector>
#include "Eigen/Dense" // 需要Eigen库支持，用来进行矩阵计算
// B-Spline基函数
double basis_function(double u, int i, int p, const std::vector<double>& knots)
{
    if (p == 0) {
        if (knots[i] <= u && u < knots[i+1]) {
            return 1.0;
        } else {
            return 0.0;
        }
    } else {
        double N1 = (u - knots[i]) / (knots[i+p] - knots[i]) * basis_function(u, i, p-1, knots);
        double N2 = (knots[i+p+1] - u) / (knots[i+p+1] - knots[i+1]) * basis_function(u, i+1, p-1, knots);
        return N1 + N2;
    }
}
// 计算NURBS曲线上某一点的坐标和一阶导数
void evaluate_single(double u, const Eigen::MatrixXd& P, const Eigen::VectorXd& w, int p, const std::vector<double>& knots, Eigen::VectorXd& r, Eigen::VectorXd& dU)
{
    int n = P.rows() - 1; // 控制点数减一
    Eigen::MatrixXd dN(p+1, n+1);
    Eigen::VectorXd N(n+1);
    // 计算基函数及其导数
    for (int i = 0; i <= n; ++i) {
        for (int j = 0; j <= p; ++j) {
            dN(j, i) = p / (knots[i+j+1] - knots[i+j]) * (basis_function(u, i+j+1, p-1, knots) - basis_function(u, i+j, p-1, knots));
        }
        N(i) = basis_function(u, i, p, knots);
    }
    // 计算点坐标
    Eigen::MatrixXd wP(n+1, 4);
    for (int i = 0; i <= n; ++i) {
        wP.row(i) << P.row(i), w(i);
    }
    r = N.transpose() * wP;
    r /= N.transpose() * w;
    // 计算一阶导数
    dU = dN * wP;
    dU /= N.transpose() * w;
    dU -= r * N.transpose() * wP * dN.transpose() * w / (N.transpose() * w * N.transpose() * w);
}
// 计算NURBS曲线弧长
double nurbs_curve_length(double a, double b, const Eigen::MatrixXd& P, const Eigen::VectorXd& w, int p, const std::vector<double>& knots, double tol)
{
    double mid = (a + b) / 2;
    Eigen::VectorXd r1, r2, dU1, dU2;
    evaluate_single(a, P, w, p, knots, r1, dU1);
    evaluate_single(b, P, w, p, knots, r2, dU2);
    double L = 0.0;
    double L1 = (r1 - r2).norm();
    double L2 = (r1 + 4 * evaluate_single(mid, P, w, p, knots, r1, dU1).tail(3) + r2).norm();
    if (std::abs(L1 - L2) < tol) {
        L = (L1 + L2) / 2;
    } else {
        L = nurbs_curve_length(a, mid, P, w, p, knots, tol/2) + nurbs_curve_length(mid, b, P, w, p, knots, tol/2);
    }
    return L;
}
int main()
{
    // 示例：计算一个NURBS曲线的弧长
    int n = 7; // 控制点数
    int p = 3; // 阶数
    std::vector<double> knots = {0, 0, 0, 0, 0.5, 1, 1, 1, 1}; // 节点矢量
    Eigen::MatrixXd P(n+1, 3); // 控制点坐标
    Eigen::VectorXd w(n+1); // 权重
    P << 0, 0, 0,
         1, 3, 0,
         2, 0, 0,
         3, 3, 0,
         4, 0, 0,
         5, 3, 0,
         6, 0, 0,
         7, 3, 0;
    w << 1, 3, 1, 3, 1, 3, 1, 3;
    double a = knots.front();
    double b = knots.back();
    double tol = 1e-6;
    double L = nurbs_curve_length(a, b, P, w, p, knots, tol);
    std::cout << "The length of the NURBS curve is: " << L << std::endl;
    return 0;
}

在 nurbs_curve_length 函数中，采用递归方法计算曲线弧长，当相邻两次积分的结果相差较小时，直接将两段的积分结果加起来，否则将这一小段再分成两半，递归调用 nurbs_curve_length 函数，直到精度满足要求为止。