构造按序经过多个点的平滑曲线，求解曲线上均匀间隔的点的坐标

绝峰·魁拔

已于 2023-04-03 22:08:33 修改

阅读量745

点赞数

文章标签：算法 python

于 2023-03-27 16:36:50 首次发布

本文链接：https://blog.csdn.net/weixin_45528232/article/details/129482835

版权

文章介绍了如何通过贝塞尔曲线构建经过多个点的平滑曲线，解释了贝塞尔曲线的特性，以及如何在Python中实现单条和多条曲线的连接。通过牛顿迭代法和高斯-勒让德公式，实现了曲线上匀速运动的计算，确保点间距离均匀分布。

摘要由CSDN通过智能技术生成

笔者首先寻找了按序经过多个点构造平滑曲线的方法，如章一；随后选择了贝塞尔曲线作为实现方法，对应的可行性分析如章二；最后展示了如何用 python 代码构造单条贝塞尔曲线、连接多条贝塞尔曲线以及实现贝塞尔曲线上的匀速运动。

一、构造平滑曲线

二、贝塞尔曲线方法构造按序经过多个点的平滑曲线的可行性

三、贝塞尔曲线及曲线上的匀速运动的 python 实现

Python 代码

一、构造平滑曲线

链接【0】：经过已知离散点画平滑曲线算法（样条曲线插值法）_画样条曲线算法_hzsjun的博客-CSDN博客

术语解释：

样条（spline curve）：在绘图术语中样条是通过一组指定点集而生成平滑曲线的柔性带。

控制点：绘制样条曲线的方法是给定一组称为控制点的坐标点，可以得到一条样条曲线。

样条曲线分为:

1 插值样条曲线(interpolate) 生成的样条曲线通过这组控制点。如牛顿插值（阶数较高时曲线不稳定）

2 逼近样条曲线(approximate) 生面的样条曲线不通过或通过部分控制点。如贝塞尔插值和 B 样条插值（笔者并不了解后一概念）

牛顿插值和平方、立方插值等只能插值函数，我们的点集可能不满足要求。

笔者最终使用贝塞尔曲线，如果对其他方法感兴趣，可以参考本章的链接【0】

二、贝塞尔曲线方法构造按序经过多个点的平滑曲线的可行性

我问了 newBing “构造按序经过多个点的平滑曲线”，可能我没有描述清楚，他给我了三次方贝塞尔曲线的答案，这是逼近样条曲线的一类解法，它不要求经过控制点。原以为无法使用贝塞尔曲线求解（因为我们要求经过所有点），幸运的是，我们找到了文章：

链接【1】：

(36条消息) 样条曲线(下)之插值问题(贝塞尔曲线、B样条和一般样条曲线插值)_b样条插值_天真的和感伤的想象家的博客-CSDN博客

由该博客，我们有：

1）贝塞尔曲线不经过中间的控制点，但必然经过两端的两个控制点

2）两个贝塞尔曲线如果平滑连接，则需要连接点与其左右相邻两端点共线

3）由2），拼接多段贝塞尔曲线，就可以让曲线经过多个指定点

因此，用贝塞尔曲线方法构造按序经过多个点的平滑曲线是可行的。

三、贝塞尔曲线及曲线上的匀速运动的 python 实现

参考：

【1】〔manim | 可视化〕你见过贝塞尔曲线吗 | 贝塞尔曲线参数方程推导_哔哩哔哩_bilibili。给定起点、终点和（零个、单个或多个）控制点，求对应贝塞尔曲线的公式

【2】曲线上的匀速运动 - 饭后温柔 - 博客园 (cnblogs.com)。博客提供了详细的公式解释，参照此博客使用高斯-勒让德公式及其对应的5阶权重表，并查询 newbing 牛顿迭代法求解函数的根后，基本解决这个问题，即实现了：能够求解某个时间质点在曲线上的位置

【3】【复合五点高斯-勒让德公式】_五点高斯勒让德求积公式_花树堆雪的博客-CSDN博客。博客提供了 5 点高斯勒让德求积公式的求积系数和求积节点

【4】匀速贝塞尔曲线运动的实现_贝塞尔曲线匀速运动_auccy的博客-CSDN博客。博客二次贝塞尔曲线匀速运动的实现并提供了 C++ 代码，笔者此时需要三次贝塞尔曲线匀速运动的 python 代码，故没有太多参考此博客，需要 C++ 实现的小伙伴可以参考。

贝塞尔曲线方程图 1 所示（图源参考【1】）：

图 1 贝塞尔曲线方程

我们可以将贝塞尔曲线写作关于 t 的函数。图 2 中将 5 段贝塞尔曲线连接起来，相同颜色代表同一贝塞尔曲线。如图 2 所示，当 t 均匀变化时，曲线上点间的距离并不是均匀变化：如蓝色段开始时，间距很小，随后间距很大。

图2 贝塞尔曲线是关于 t 的函数, 其中 t 的范围为 [0,1]，本示例 t 均匀增长

我们想要得到贝塞尔曲线上点间距离分布均匀的离散点，我们参考了【2】。贝塞尔曲线的两个端点，可以分别称之为起点和终点，【2】中的思路是：

1）利用贝塞尔曲线函数的导数求积分，得到从起点沿着曲线到曲线上某一点的距离的、关于 t 的函数 L(t)。这里的积分可以用【3】中的牛顿勒让德求积近似公式完成

2）然后给定任意距离 s，用牛顿迭代法可以求解 L(t) = s 等式中的 t，将 t 代回贝塞尔曲线，就可以得到贝塞尔曲线上距离起点距离为 s 的点的坐标。

3）给定间隔均匀的距离列表（如 [0.3, 0.6, 0.9] 间隔均为 0.3），依次用 2）中的步骤求解对应的坐标，就可以得到贝塞尔曲线上距离均匀的点。如图 3 和图 4（图 3 中，仅单段贝塞尔曲线上距离均匀，图 4 中，多段贝塞尔曲线上距离均匀）

由于求积分时使用了牛顿勒让德近似公式，因此得到的关于 t 的距离函数 L(t) 并不完全精确，所以按1）~ 3）的方法求解得到的点也不是完全距离均匀的。

图3 单个贝塞尔曲线上距离相同的点（不同颜色表示不同的贝塞尔曲线）

图 3 中，同一贝塞尔曲线上的中间部分的点间距相同，但两端的点由于曲线长度并不能整除设定的间距，所以间距略有不同。在代码上略修改，可以得到多段连接的贝塞尔曲线上距离相同的点，如图4.

图4 相连贝塞尔曲线上距离相同的点，黄色点为不同贝塞尔曲线的连接点

Python 代码

文件：toys/bizierAndMore.py at main · TuZhengzhou/toys (github.com)

以下代码与文件中内容一致

from scipy.special import comb
import numpy as np
import matplotlib.pyplot as plt
from typing import List
# from tools import 

def QT(controlPoints, t):
    """
    返回贝塞尔曲线上 t 对应的点的坐标
    """
    n = len(controlPoints)-1
    Bt = np.zeros(2, np.float64)
    for i in range(len(controlPoints)): # 四个点使用 3 的组合, 即 C(3,0), C(3,1), C(3,2), C(3,3)
        Bt = Bt + comb(n,i) * np.power(1-t,n-i) * np.power(t,i) * np.array(controlPoints[i])
    return Bt

def QT_(controlPoints, t):
    """
    Q(t) 函数的一阶导数
    """
    n = len(controlPoints)-1
    Bt = np.zeros(2, np.float64)
    for i in range(len(controlPoints)): 
        # Bt = Bt + comb(n,i) * np.power(1-t,n-i) * np.power(t,i) * np.array(controlPoints[i])
        c1 = - (n-i) * np.power(1-t,n-i-1) * np.power(t,i)
        c2 = i * np.power(1-t,n-i) * np.power(t,i-1)
        Bt = Bt + comb(n,i) * (c1 + c2) * np.array(controlPoints[i])
    return Bt

"""高斯勒让德积分公式的求积系数和求积节点"""
W = [0.5688888888888889, 0.4786286704993665, 0.4786286704993665, 0.2369268850561891, 0.2369268850561891]
X = [0.0000000000000000, -0.5384693101056831, 0.5384693101056831, -0.9061798459386640, 0.9061798459386640]

def LT(controlPoints, t) -> float:
    """
    返回 点 QT(t) 到贝塞尔曲线起点的距离
    """
    ret = float(0)
    n = W.__len__()
    t_half = t/2
    for i in range(n):
        qt_ = QT_(controlPoints, t_half*(X[i]+1))
        ret += W[i] * np.sqrt(qt_[0]*qt_[0] + qt_[1]*qt_[1])
    ret *= t_half
    return ret

def LT_(controlPoints, t):
    """
    LT() 的一阶导数
    """
    qt_ = QT_(controlPoints, t)
    return np.sqrt(qt_[0]*qt_[0] + qt_[1]*qt_[1])

def newton(controlPoints, start_x, target_y, tol) -> float:
    """
    牛顿迭代法求解满足 target_y = LT(t) 的 t
    start_x: t 的初始值
    target_y: 此处指目标距离
    tol: 与 target_y 的最大允许误差
    """
    fx = LT(controlPoints, start_x)
    # print("start: ", start_x)
    while abs( fx - target_y ) > tol:
        dfx = LT_(controlPoints, start_x)
        start_x = start_x - (fx - target_y) / dfx
        fx = LT(controlPoints, start_x)
        # print(start_x)
    return start_x

def getInterpolationPoints(controlPoints, tList):
    """
    给出 t 的列表，返回贝塞尔曲线上对应点的坐标
    """
    n = len(controlPoints)-1
    interPoints = []
    for t in tList:
        Bt = np.zeros(2, np.float64)
        for i in range(len(controlPoints)):
            Bt = Bt + comb(n,i) * np.power(1-t,n-i) * np.power(t,i) * np.array(controlPoints[i])
        # Bt = [t*50, LT(controlPoints, t)]
        # Bt = [t*50, LT_(controlPoints, t)]
        interPoints.append(list(Bt))

    return interPoints


def getControlPointList(pointsArray):
    """
    求解贝塞尔曲线的控制点
    这个算法可以自行设计, 控制点会决定曲线的形状
    """
    points = np.array(pointsArray)
    index = points.shape[0] - 2

    res = []
    for i in range(index):

        tmp = points[i:i+3]
        p1 = tmp[0]
        p2 = tmp[1]
        p3 = tmp[2]
        
        l12 = np.sqrt(np.sum((p2 - p1)**2))
        l23 = np.sqrt(np.sum((p3 - p2)**2))
        p12_norm = (p2 - p1) / l12
        p23_norm = (p3 - p2) / l23
        e = (p12_norm + p23_norm) / np.sqrt(np.sum((p12_norm + p23_norm)**2))
        c1 = p2 - e * l12 * 0.2
        c2 = p2 + e * l23 * 0.2
        res.append(c1)
        res.append(c2)

    pFirst = points[0] + 0.1*(res[0] - points[0])
    pEnd = points[-1] + 0.1*(res[-1] - points[-1])
    res.insert(0,pFirst)
    res.append(pEnd)

    return np.array(res)

def getPositionsFromDist(dists: List[float], points: List[List[float]], controlP: np.ndarray) -> List[List[float]]:
    """
    我们按照一个划定的路线行进，本函数输入在路线上行进的距离，返回与距离一一对应的坐标
    
    dists 表示距离列表
    points 和 controlP 用于限制路线
    
    返回坐标列表
    """
    l = len(points) - 1
    figure = plt.figure()
    j = 0
    extraLen = 0
    retPoints = []
    for i in range(l):
        controlPoints = np.array([points[i], controlP[2*i], controlP[2*i+1], points[i+1]])
        total_len = LT(controlPoints, 1)
        
        tList = []
        for dist in dists[j:]:
            if dist - extraLen > total_len:
                j += len(tList)
                extraLen += total_len
         
                # plt.scatter(range(len(tList)),tList)
                retPoints += getInterpolationPoints(controlPoints, tList)
                break
            else:
                target_y = dist - extraLen
                start_t = target_y / total_len
                t = newton(controlPoints, start_t, target_y, 0.0000000001)
                tList.append(t)
        
    x = np.array(retPoints)[:,0]
    y = np.array(retPoints)[:,1]
    plt.scatter(x,y)
    plt.plot(x, y)
    plt.scatter(np.array(points)[:,0],np.array(points)[:,1])
    plt.show()
    return retPoints
        
def example1():
    """
        贝塞尔曲线是关于 t 的函数, 其中 t 的范围为 [0,1]
        本示例 t 均匀增长, 画出 t 均匀增长时曲线上点的位置
    """
    points = [[1,1],[3,6],[6,3],[8,0],[11,6],[12,12]]
    controlP = getControlPointList(points)
    l = len(points) - 1
    figure = plt.figure()
    t = np.linspace(0,1,12)
    for i in range(l):
        p = np.array([points[i], controlP[2*i], controlP[2*i+1], points[i+1]])
        interPoints = getInterpolationPoints(p, t)
        x = np.array(interPoints)[:,0]
        y = np.array(interPoints)[:,1]
        plt.scatter(x,y)
        plt.plot(x, y)
            
    plt.scatter(np.array(points)[:,0],np.array(points)[:,1],color='gray')
    plt.show()
    
def example2():
    """
        贝塞尔曲线是关于 t 的函数, 其中 t 的范围为 [0,1]
        
        本示例首先求出贝塞尔曲线的长度 s 随 t 变化的公式 LT(controlPoints, t), 当 t 为 0 时该函数值为 0, 当 t 为 1 时该函数值为贝塞尔曲线全长
        
        我们用牛顿迭代法求解与指定 s 对应的 t, 再使用求解出的 t 算出坐标, 从而求解与指定长度 s 对应的曲线上的坐标
    """
    points = [[1,1],[3,6],[6,3],[8,0],[11,6],[12,12]]
    controlP = getControlPointList(points)
    l = len(points) - 1
    
    figure = plt.figure()
    for i in range(l):
        p = np.array([points[i], controlP[2*i], controlP[2*i+1], points[i+1]])
        
        total_len = LT(p, 1)
        target_ys = np.linspace(0, total_len, max(2, math.ceil(total_len/0.5)))
        print(target_ys)
        
        ts = []
        for target_y in target_ys:
            start_x = target_y / total_len
            x = newton(p, start_x, target_y, 0.0000001)
            ts.append(x)

        interPoints = getInterpolationPoints(p, ts)
        x = np.array(interPoints)[:,0]
        y = np.array(interPoints)[:,1]
        plt.scatter(x,y)
        plt.plot(x,y)
            
    plt.scatter(np.array(points)[:,0],np.array(points)[:,1],color='gray')
    plt.show()

def example3():
    """
    是 example2 的扩展版, 我们连接了多段 贝塞尔曲线, 这意味着我们给出的长度是由多段贝塞尔曲线相加得到的（而不是在单段贝塞尔曲线上）
    """
    dists = np.linspace(0,30,68)
    print(dists)
    points = [[1,1],[3,6],[6,3],[8,0],[11,6],[12,12]]
    controlP = getControlPointList(points)
    retPoints = getPositionsFromDist(dists, points, controlP)
    print(retPoints.__len__())

import math
if __name__ == '__main__':
    example1()      # 对应 图 2
    example2()      # 对应 图 3
    example3()      # 对应 图 4

绝峰·魁拔

关注

0
点赞
踩
9

收藏

觉得还不错? 一键收藏
0
评论
构造按序经过多个点的平滑曲线，求解曲线上均匀间隔的点的坐标

笔者首先寻找了按序经过多个点构造平滑曲线的方法，如章一；随后选择了贝塞尔曲线作为实现方法，对应的可行性分析如章二；最后展示了如何用 python 代码构造单条贝塞尔曲线、连接多条贝塞尔曲线以及实现贝塞尔曲线上的匀速运动。
复制链接

扫一扫