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

Simpson 积分方法是一种数值积分方法，可以用于计算曲线的弧长。它的基本思想是将曲线分成若干小段，对每一小段采用 Simpson 公式进行数值积分，然后将各小段的积分结果加起来得到整个曲线的积分结果。
具体地，设曲线 C C C 由参数方程 r ( t ) = ( x ( t ) , y ( t ) ) \mathbf{r}(t)=(x(t), y(t)) r(t)=(x(t),y(t)) 给出， t ∈ [ a , b ] t \in [a, b] t∈[a,b]， n n n 是分段数，则将 t t t 分成 n n n 段，即 t 0 = a , t 1 , … , t n = b t_0=a, t_1, \ldots, t_n=b t0=a,t1,…,tn=b，每段长度为 h = ( b − a ) / n h=(b-a)/n h=(b−a)/n，则 Simpson 公式可以表示为：
∫ t i t i + 2 ∣ r ′ ( t ) ∣ d t ≈ h 3 ( ∣ r ′ ( t i ) ∣ + 4 ∣ r ′ ( t i + 1 ) ∣ + ∣ r ′ ( t i + 2 ) ∣ ) ( i = 0 , 2 , 4 , … , n − 2 ) \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) ∫titi+2∣r′(t)∣dt≈3h(∣r′(ti)∣+4∣r′(ti+1)∣+∣r′(ti+2)∣)(i=0,2,4,…,n−2)
这个公式是对 [ t i , t i + 2 ] [t_i, t_{i+2}] [ti,ti+2] 区间上的函数 ∣ r ′ ( t ) ∣ |\mathbf{r}'(t)| ∣r′(t)∣ 进行数值积分的，近似结果可以看作是这个区间上函数曲线下面的面积，即弧长的近似值。将所有区间的积分结果加起来即得到整个曲线的积分结果。
对于 NURBS 曲线，我们可以计算其导数 r ′ ( t ) \mathbf{r}'(t) 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 函数，直到精度满足要求为止。

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

matlab蒙特卡罗方法求体积_实验二-蒙特卡罗方法计算三维体积
班级: 学号: 姓名: 实验时间: 2014 年月日实验项目实验二蒙特卡罗方法计算三维体积所属课程数学实验实验目的了解蒙特卡罗方法的原理,掌握随机数使用技术. 实验内 ...
利用NURBS曲线进行点云曲面拟合算法
文章目录介绍 NURBS曲线 C++实现思路代码实现读取点云数据对点云进行预处理创建曲面模型将曲面模型转换为NURBS曲面完整代码 opennurbs.h说明 vs2019安装OpenN ...
【四足机器人--摆动相足端位置速度轨迹规划】（4.1）FootSwingTrajectory（bezier曲线计算脚的摆动轨迹）代码解析
系列文章目录提示:这里可以添加系列文章的所有文章的目录,目录需要自己手动添加 TODO:写完再整理文章目录系列文章目录前言一.FootSwingTrajectory(bezier曲线)的内容 ...
Algorithm：机械优化设计的数学模型简介、常用优化方法、优化计算工具简介之详细攻略
Algorithm:机械优化设计的数学模型简介.常用优化方法.优化计算工具简介之详细攻略目录机械设计中基于算法模型的机械优化设计 1.优化设计的数学模型
计算nurbs曲率NURBS曲线的曲率计算
计算nurbs曲率NURBS曲线的曲率计算 https://download.csdn.net/download/yishang44/10383324 NURBS曲线的计算和求导 ,主要看的书籍是&l ...
NURBS曲线的曲率计算
NURBS曲线的曲率计算这两天做课题需要用到NURBS曲线,仔细探究了NURBS曲线的计算和求导等,主要看的书籍是<The NURBS book>,代码主要参考nurbs工具箱-M语言, ...
通俗易懂的MonteCarlo积分方法(七)
通俗易懂的MonteCarlo积分方法(七) 一.设计的基本GUIGUIGUI格式menteCarlo.figmenteCarlo.figmenteCarlo.fig 二.设计MatlabMatlab ...
运用数学软件matlab求无穷积分,matlab积分的计算及其简单应用论文.doc
积分的计算及其简单应用摘要:本文简要的概述了MATLAB 在高等数学中积分的计算及应用:利用MATLAB 中符号积分和数值积分的命令,计算定积分和不定积分.同时,也可以通过这些命令来解决一些实际问题 ...
第一型曲线积分与第一型曲面积分、第二型曲线积分与格林公式
提示:本文的适用对象为已修过<微积分A1>的非数学系学生,文中题型方法为个人总结,为个人复习使用.部分理解虽然不太严谨,但对于解题的实用性较强.若有疏漏or错误,欢迎批评指正. 一.关于第 ...

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

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

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

最新文章

热门文章