辛普森法则（Simposon rule）求参数曲线弧长 Python （NURBS Butterfly 为例）

NURBS

非均匀有理样条曲线 （由控制点Pi， 权重wi，节点向量u，阶数k定义）

辛普森法则

已知参数曲线某段弧长的两个端点[a, b]，求这段弧的弧长


在下面的python代码中c=(a+b)/2

Python 代码实现

from geomdl import linalg, NURBS, compatibility
import numpy as np
import math

# 求得曲线一阶导
def first_derivative(a, crv):
    ders = crv.derivatives(a, order=1)
    der1 = ders[1]
    A = math.sqrt(der1[0] ** 2 + der1[1] ** 2 + der1[2] ** 2)
    return A


def simpson(a, b, crv):
    c = a + (b - a) / 2
    s = (first_derivative(a, crv) + 4 * first_derivative(c, crv) + first_derivative(b, crv)) * (b - a) / 6
    return s


def arc_length(a, b, crv, A, epsilon=0.0000001):
    c = a + (b - a) / 2
    L = simpson(a, c, crv)
    R = simpson(c, b, crv)

    if abs(L + R - A) <= 10 * epsilon:
        length = L + R + (L + R - A) / 10
    else:
        length = arc_length(a, c, crv, epsilon / 2, L) + arc_length(c, b, crv, epsilon / 2, R)
    return length


def arc_length1(us, ue, crv):
    LE = 0
    for i in range(0, 100000 - 1):
        uk = us + i * (ue - us) / 100000
        uk1 = us + (i + 1) * (ue - us) / 100000
        der = crv.evaluate_single(uk)
        der1 = crv.evaluate_single(uk1)
        Length = math.sqrt((der[0] - der1[0]) ** 2 + (der[1] - der1[1]) ** 2 + (der[2] - der1[2]) ** 2)
        LE += Length
    return LE


if __name__ == '__main__':
	# butterfly NURBS曲线参数
	crv = NURBS.Curve()
    crv.degree = 3
    ctrlpts = [[0, 54.493, 52.139], [0, 55.507, 52.139], [0, 56.082, 49.615], [0, 56.780, 44.971], [0, 69.575, 51.358],
               [0, 77.786, 58.573], [0, 90.526, 67.081], [0, 105.973, 63.801], [0, 100.400, 47.326],
               [0, 94.567, 39.913],
               [0, 92.369, 30.485], [0, 83.440, 33.757], [0, 91.892, 28.509], [0, 89.444, 20.393], [0, 83.218, 15.446],
               [0, 87.621, 4.830], [0, 80.945, 9.267], [0, 79.834, 14.535], [0, 76.074, 8.522], [0, 70.183, 12.550],
               [0, 64.171, 16.865], [0, 59.993, 22.122], [0, 55.680, 36.359], [0, 56.925, 24.995], [0, 59.765, 19.828],
               [0, 54.493, 14.940], [0, 49.220, 19.828], [0, 52.060, 24.994], [0, 53.305, 36.359], [0, 48.992, 22.122],
               [0, 44.814, 16.865], [0, 38.802, 12.551], [0, 32.911, 8.521], [0, 29.152, 14.535], [0, 28.040, 9.267],
               [0, 21.364, 4.830], [0, 25.768, 15.447], [0, 19.539, 20.391], [0, 17.097, 28.512], [0, 25.537, 33.750],
               [0, 16.602, 30.496], [0, 14.199, 39.803], [0, 8.668, 47.408], [0, 3.000, 63.794], [0, 18.465, 67.084],
               [0, 31.197, 58.572], [0, 39.411, 51.358], [0, 52.204, 44.971], [0, 52.904, 49.614], [0, 53.478, 52.139],
               [0, 54.493, 52.139]]
	w = [1.00, 1.00, 1.00, 1.20, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 2.00, 1.00, 1.00, 5.00, 3.00, 1.00,
         1.10, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.000, 1.00, 1.10,
         1.00, 3.00, 5.00, 1.00, 1.00, 2.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.20, 1.00, 1.00, 1.00]
    crv.ctrlptsw = compatibility.combine_ctrlpts_weights(ctrlpts, w)
    crv.knotvector = [0, 0, 0, 0, 0.0083, 0.0150, 0.0361, 0.0855, 0.1293, 0.1509, 0.1931, 0.2273, 0.2435, 0.2561,
                      0.2692, 0.2889, 0.3170, 0.3316, 0.3482, 0.3553, 0.3649, 0.3837, 0.4005, 0.4269, 0.4510,
                      0.4660, 0.4891, 0.5000, 0.5109, 0.5340, 0.5489, 0.5731, 0.5994, 0.6163, 0.6351, 0.6447,
                      0.6518, 0.6683, 0.6830, 0.7111, 0.7307, 0.7439, 0.7565, 0.7729, 0.8069, 0.8491, 0.8707,
                      0.9145, 0.9639, 0.9850, 0.9917, 1, 1, 1, 1] 
a = 0.01
b = 0.1
A = simpson(a, b, crv)
Length =  arc_length(a, b, crv, A, epsilon=0.0000001)
Length1 = arc_length1(a, b, crv)
print(Length)
print(Length1)  

print
36.70465969584353
33.49012043159443

不知道为什么两种方法求出来的弧长不一样，各位如果知道为什么，请在评论区或私信告诉我

geomdl 样条曲线库

geomdl 样条曲线库太好用了，具体见
https://nurbs-python.readthedocs.io/en/5.x/index.html