路径平滑优化算法--B样条(B-spline)路径平滑算法
路径平滑优化算法–B样条(B-spline)路径平滑算法
文章目录
- 路径平滑优化算法--B样条(B-spline)路径平滑算法
- 1.贝塞尔曲线的局限性
- 2.B样条曲线的提出
- 3. 关键定义
- 3.1 节点向量的具体示例
- 4.曲线方程
- 5.算法步骤
- 步骤1:准备输入数据和参数
- 步骤2:参数化路径点
- 步骤3:构建基矩阵并求解控制点
- 步骤4:生成平滑曲线
- 6. B样条路径平滑计算实例(逐步详解)
- 6.1 问题描述
- 步骤 1:确定参数
- 步骤 2:对路径点进行参数化(Chord-length)
- 步骤 3:计算 B 样条基函数矩阵 AAA
- 步骤 4:最小二乘求解控制点 P\mathbf{P}P
- **step 0 准备数据**
- **step 1**:计算基函数矩阵 AAA
- **step 2:解线性系统 A⋅Px=QxA \cdot P_x = Q_xA⋅Px=Qx**
- **step 3:数值求解(用代码计算)**
- 步骤 5:构造曲线并采样
- 7.为什么 B 样条的“构造 + 采样” = 路径平滑?
- 7.1 构造形式的回顾
- 7.2 为什么这样可以实现平滑路径优化?
- 7.2.1. 基函数提供“天然平滑基底”
- 7.2.2. 控制点的解来自于“最小二乘优化问题”
- 7.2.3. 节点向量控制曲线的局部性和形状
- 7.3 将这套方法理解为“正则化的逼近问题”
- 8.python实现代码
在机器人学、计算机图形和自主导航领域,从离散路径点生成平滑连续路径至关重要,以确保高效安全的运动。B样条路径平滑算法通过将锯齿状或有角路径转换为流畅曲线来解决这一需求。该方法利用B样条曲线的数学特性,减少突然变化,从而降低机器人或车辆系统的机械应力。
本文以博客形式详细一步一步解释B样条路径平滑算法。从基础概念开始,到算法原理,然后分解成可操作步骤。读者将清晰理解如何实现和应用此技术。假设读者对线性代数和参数曲线有基本了解,但会提供清晰解释。
1.贝塞尔曲线的局限性
尽管贝塞尔曲线以其简单性和良好的整体连续性被广泛应用于曲线建模,但在实际应用中仍存在以下关键缺陷:
- 阶数与控制点数量强绑定
贝塞尔曲线的阶数由控制点数量决定,若有 n+1n+1n+1 个控制点,则曲线为 nnn 次多项式。例如,若使用 100 个控制点,生成的为 99 次多项式曲线,其导数将具有多达 98 个极值点,极易引发 龙格现象(Runge Phenomenon),导致曲线剧烈震荡、不可控。 - 拼接复杂性高
多段贝塞尔曲线的拼接需显式约束控制点的位置以维持高阶连续性(如 C1,C2C^1, C^2C1,C2),增加建模难度,尤其在自动生成路径的场景中效果不理想。 - 缺乏局部可调性
贝塞尔曲线的另一个重要缺陷在于:全局控制性强但局部修改困难。调整任意一个控制点会影响整个曲线形状,缺乏“局部性”限制,不利于局部优化或编辑。
2.B样条曲线的提出
为克服上述问题,Gordon 和 Riesenfeld 等人在 1972 年提出了 B样条曲线(Basis Spline Curve) 的概念。B样条保留了贝塞尔曲线的优点(如控制点可视性、光滑性),同时具有以下优势:
- 低阶分段表达,抑制震荡
B样条曲线以一组低次多项式(通常为三次)分段定义,避免使用高阶全局多项式,从而有效抑制震荡行为。 - 自然拼接,高阶连续
相邻曲段之间可自动保证一定阶数的连续性(如三次B样条保证 C2C^2C2 连续),无需手动调整。 - 局部可修改性
每个控制点只影响其支持范围内的曲线段,实现“牵一发而动局部”的效果,大大提升了建模灵活性和效率。
上图直观对比了 Bézier 曲线与 B 样条曲线在控制点局部修改时的响应差异:
左图(Bézier 曲线):修改中间一个控制点后,整条曲线发生了全局变化,体现出“牵一发而动全身”的特性。
右图(B 样条曲线):相同修改仅影响部分曲段,其他部分保持不变,体现出良好的局部可控性。
3. 关键定义
-
控制点:一组点 {Pi}\{\mathbf{P}_i\}{Pi}(i从0到ni 从 0 到 ni从0到n),影响曲线形状但不一定位于曲线上。
-
度数 (k):曲线段的多项式度。通常 k=3k=3k=3 用于三次B样条,确保二阶 (G2G^2G2) 连续性。
-
节点向量:非递减序列 t=[t0,t1,…,tm]\mathbf{t} = [t_0, t_1, \dots, t_{m}]t=[t0,t1,…,tm],其中 m=n+k+1m = n + k + 1m=n+k+1。节点定义参数区间。
-
基函数:B样条基函数 Ni,k(u)N_{i,k}(u)Ni,k(u) 通过Cox-de Boor公式递归计算:
- Ni,0(u)={1如果 ti≤u<ti+10否则N_{i,0}(u) = \begin{cases} 1 & \text{如果 } t_i \leq u < t_{i+1} \\ 0 & \text{否则} \end{cases}Ni,0(u)={10如果 ti≤u<ti+1否则
- Ni,k(u)=u−titi+k−tiNi,k−1(u)+ti+k+1−uti+k+1−ti+1Ni+1,k−1(u)N_{i,k}(u) = \frac{u - t_i}{t_{i+k} - t_i} N_{i,k-1}(u) + \frac{t_{i+k+1} - u}{t_{i+k+1} - t_{i+1}} N_{i+1,k-1}(u)Ni,k(u)=ti+k−tiu−tiNi,k−1(u)+ti+k+1−ti+1ti+k+1−uNi+1,k−1(u)
- 这些函数确保曲线是分段多项式,并满足单位划分性质: ∑iNi,k(u)=1\sum_i N_{i,k}(u) = 1∑iNi,k(u)=1
3.1 节点向量的具体示例
我们希望构造一个“开放均匀节点向量(Open Uniform Knot Vector)”,使得:
-
曲线在首尾“逼近”端点(但不一定插值)
-
控制曲线的分段与过渡区
-
基函数计算具备数值稳定性
-
nnn:控制点数量 − 1,即曲线的“阶数减1”
-
kkk:B样条曲线的阶数(比如三次B样条 k=3k=3k=3)
-
节点向量 t={t0,t1,…,tn+k+1}t = \{t_0, t_1, \dots, t_{n+k+1} \}t={t0,t1,…,tn+k+1},共 n+k+2n+k+2n+k+2 个值
-
示例:
-
控制点个数:5 个 ⇒ n=4n = 4n=4
-
阶数:3(即三次B样条)
-
所需节点向量长度:
长度=n+k+2=4+3+2=9\text{长度} = n + k + 2 = 4 + 3 + 2 = 9长度=n+k+2=4+3+2=9
构造方式(Open Uniform)
构造规则如下:
- 开头重复 k+1=4k+1 = 4k+1=4个最小值(如0)
- 结尾重复 k+1=4k+1 = 4k+1=4个最大值(如1)
- 中间插入均匀分布的内部节点(如0.5)
t=[0,0,0,0]⏟起点重复 k+1次,中间节点=[0.5],[1,1,1,1]⏟终点重复 k+1次t = \underbrace{[0, 0, 0, 0]}_{\text{起点重复 }k+1\text{ 次}}, \quad \text{中间节点} = [0.5], \quad \underbrace{[1, 1, 1, 1]}_{\text{终点重复 }k+1\text{ 次}}t=起点重复 k+1 次[0,0,0,0],中间节点=[0.5],终点重复 k+1 次[1,1,1,1]
最终结果:
t=[0,0,0,0,0.5,1,1,1,1]t = [0, 0, 0, 0, 0.5, 1, 1, 1, 1]t=[0,0,0,0,0.5,1,1,1,1]
📝 说明:
- 起始段:保证曲线从第一个控制点开始
- 结尾段:保证曲线向最后一个控制点靠拢
- 中间节点数量=n−k+1=2中间节点数量 = n−k+1=2中间节点数量=n−k+1=2,其中插入1个内部节点
图片展示了Cubic B-Spline Basis Functions(三次B样条基函数)在开放均匀节点向量条件下的形状变化:
每条曲线表示一个基函数 𝑁𝑖,3(𝑢),仅在局部范围非零;各基函数平滑连接,整体构成分段多项式族;用于组合控制点生成光滑曲线。
4.曲线方程
B样条曲线 C(u)\mathbf{C}(u)C(u) 表示为:
C(u)=∑i=0nNi,k(u)Pi,u∈[tk,tn+1]\mathbf{C}(u) = \sum_{i=0}^{n} N_{i,k}(u) \mathbf{P}_i, \quad u \in [t_k, t_{n+1}]C(u)=∑i=0nNi,k(u)Pi,u∈[tk,tn+1]
导数如速度 C′(u)\mathbf{C}'(u)C′(u) 和加速度 C′′(u)\mathbf{C}''(u)C′′(u) 可类似求得,用于计算曲率:
κ(u)=∣C′(u)×C′′(u)∣∣C′(u)∣3\kappa(u) = \frac{|\mathbf{C}'(u) \times \mathbf{C}''(u)|}{|\mathbf{C}'(u)|^3}κ(u)=∣C′(u)∣3∣C′(u)×C′′(u)∣
这些特性使B样条适合路径平滑,因为允许局部调整而不影响整个曲线。
5.算法步骤
算法将离散路径点{Qj}\{\mathbf{Q}_j\}{Qj}(如A*算法生成的)逼近为平滑B样条曲线。目标是最小化与原路径偏差,同时确保平滑,通常通过最小二乘拟合实现。该方法在需要避障和最小曲率变化的环境中特别有效。
变体包括均匀B样条(简单)或非均匀有理B样条 (NURBS)(加权控制和精确圆锥截面)。过程强调逼近而非插值,以避免尖锐转弯处的振荡。
这里以三次B样条 (k=3k=3k=3) 为例,针对2D路径,结构化顺序说明。每步包括数学依据和实际考虑。
步骤1:准备输入数据和参数
- 输入:离散路径点 ${\mathbf{Q}j = (x_j, y_j)}{j=0}^m ,,,m+1$ 个点。
- 参数选择:
- 选择度数 k=3k=3k=3 以确保G2G^2G2连续性。
- 选择控制点数 n+1n+1n+1,通常 n≈m/2n ≈ m/2n≈m/2 以平衡保真度和平滑性。
- 定义节点向量。对于准均匀节点:端点重复k次(例如 t0=t1=...=tk−1=0,内部均匀分布,tn+2=...=tm=1t_0 = t_1 = ... = t_{k-1} = 0,内部均匀分布,t_{n+2} = ... = t_m = 1t0=t1=...=tk−1=0,内部均匀分布,tn+2=...=tm=1)。
- 依据:适当参数化防止数值不稳定。端点重复确保曲线若需可起始于首末控制点。
- 实际提示:使用Python的SciPy库生成节点,避免手动错误。
步骤2:参数化路径点
- 为每个 Qj\mathbf{Q}_jQj 分配参数值 uju_juj。
- 方法:常用弦长参数化:
- 计算累积距离:d0=0,dj=dj−1+∣Q∗j−Q∗j−1∣d_0 = 0,d_j = d_{j-1} + |\mathbf{Q}*j - \mathbf{Q}*{j-1}|d0=0,dj=dj−1+∣Q∗j−Q∗j−1∣。
- 归一化:uj=dj/dmu_j = d_j / d_muj=dj/dm,得到 uj∈[0,1]u_j ∈ [0, 1]uj∈[0,1]。
- 依据:此方法反映几何间距,比均匀参数化提高拟合精度,后者可能扭曲不均匀点。
- 实际提示:对于高度不规则路径,考虑向心参数化 (uj∝√dj)(u_j ∝ √d_j)(uj∝√dj) 以减少过冲。
步骤3:构建基矩阵并求解控制点
- 形成方程系统:Qj≈∑i=0nNi,k(uj)Pi\mathbf{Q}_j \approx \sum_{i=0}^n N_{i,k}(u_j) \mathbf{P}_iQj≈∑i=0nNi,k(uj)Pi。
- 矩阵形式:创建 (m+1) × (n+1) 矩阵 A,其中 Aj,i=Ni,k(uj)A_{j,i} = N_{i,k}(u_j)Aj,i=Ni,k(uj)。
- 求解:使用最小二乘: min∥AP−Q∥2⟹P=(ATA)−1ATQ\min \|\mathbf{A} \mathbf{P} - \mathbf{Q}\|^2 \implies \mathbf{P} = (\mathbf{A}^T \mathbf{A})^{-1} \mathbf{A}^T \mathbf{Q}min∥AP−Q∥2⟹P=(ATA)−1ATQ
- 依据:逼近允许平滑偏差。伪逆处理超定系统 (m>nm > nm>n)。
- 实际提示:代码中使用 numpy.linalg.lstsq 高效求解。若 ATAA^T AATA病态,添加正则化(如岭回归)。
步骤4:生成平滑曲线
- 在 [0,1][0, 1][0,1] 上均匀采样 usu_sus(例如1000点)。
- 计算 C(us)=∑i=0nNi,k(us)Pi\mathbf{C}(u_s) = \sum_{i=0}^n N_{i,k}(u_s) \mathbf{P}_iC(us)=∑i=0nNi,k(us)Pi。
- 依据:采样将连续曲线离散化,用于实际如运动控制器。
- 实际提示:实时应用中预计算基函数以提高效率。
6. B样条路径平滑计算实例(逐步详解)
6.1 问题描述
给定一组二维路径点:
Q={Q0=(0,0),Q1=(1,1),Q2=(2,0),Q3=(3,1),Q4=(4,2)}\mathbf{Q} = \left\{ \begin{aligned} &Q_0 = (0, 0), \\ &Q_1 = (1, 1), \\ &Q_2 = (2, 0), \\ &Q_3 = (3, 1), \\ &Q_4 = (4, 2) \end{aligned} \right\}Q=⎩⎨⎧Q0=(0,0),Q1=(1,1),Q2=(2,0),Q3=(3,1),Q4=(4,2)⎭⎬⎫
目标:用三次 B 样条平滑拟合该路径,生成一条连续可导的曲线,并逐步计算每一阶段的具体值。
步骤 1:确定参数
曲线阶数
k=3(三次B样条)k = 3 \quad \text{(三次B样条)}k=3(三次B样条)
控制点个数
n+1=5(即与路径点个数相同)n + 1 = 5 \quad (\text{即与路径点个数相同})n+1=5(即与路径点个数相同)
节点向量(开放均匀)
t={0,0,0,0,0.5,1,1,1,1}t = \{0, 0, 0, 0, 0.5, 1, 1, 1, 1\}t={0,0,0,0,0.5,1,1,1,1}
- 前后各重复k=3k=3k=3 次端点值
- 中间设置一个内部节点 t4=0.5t_4 = 0.5t4=0.5
步骤 2:对路径点进行参数化(Chord-length)
计算路径点间距离
d1=∥Q1−Q0∥=(1−0)2+(1−0)2=2≈1.4142d2=∥Q2−Q1∥=(2−1)2+(0−1)2=2≈1.4142d3=∥Q3−Q2∥=(3−2)2+(1−0)2=2≈1.4142d4=∥Q4−Q3∥=(4−3)2+(2−1)2=2≈1.4142\begin{aligned} &d_1 = \|Q_1 - Q_0\| = \sqrt{(1-0)^2 + (1-0)^2} = \sqrt{2} \approx 1.4142 \\ &d_2 = \|Q_2 - Q_1\| = \sqrt{(2-1)^2 + (0-1)^2} = \sqrt{2} \approx 1.4142 \\ &d_3 = \|Q_3 - Q_2\| = \sqrt{(3-2)^2 + (1-0)^2} = \sqrt{2} \approx 1.4142 \\ &d_4 = \|Q_4 - Q_3\| = \sqrt{(4-3)^2 + (2-1)^2} = \sqrt{2} \approx 1.4142 \end{aligned}d1=∥Q1−Q0∥=(1−0)2+(1−0)2=2≈1.4142d2=∥Q2−Q1∥=(2−1)2+(0−1)2=2≈1.4142d3=∥Q3−Q2∥=(3−2)2+(1−0)2=2≈1.4142d4=∥Q4−Q3∥=(4−3)2+(2−1)2=2≈1.4142
总长度=d1+d2+d3+d4=4×2≈4×1.4142=5.6568总长度=d_1+d_2+d_3+d_4=4×2≈4×1.4142=5.6568总长度=d1+d2+d3+d4=4×2≈4×1.4142=5.6568
再根据以下公式归一化路径点参数:
uj=累积距离到第j个点5.6568u_j = \frac{\text{累积距离到第}j\text{个点}}{5.6568}uj=5.6568累积距离到第j个点
结果:
累加距离并归一化为参数 uju_juj
D0=0D1=1.4142D2=2.8284D3=4.2426D4=5.6568⇒u0=0/5.6568=0.0000u1=1.4142/5.6568=0.2500u2=2.8284/5.6568=0.5000u3=4.2426/5.6568=0.7500u4=1.0000\begin{aligned} &D_0 = 0 \\ &D_1 = 1.4142 \\ &D_2 = 2.8284 \\ &D_3 = 4.2426 \\ &D_4 = 5.6568 \end{aligned} \quad\Rightarrow\quad \begin{aligned} &u_0 = 0 / 5.6568 = 0.0000 \\ &u_1 = 1.4142 / 5.6568 = 0.2500 \\ &u_2 = 2.8284 / 5.6568 = 0.5000 \\ &u_3 = 4.2426 / 5.6568 = 0.7500 \\ &u_4 = 1.0000 \end{aligned}D0=0D1=1.4142D2=2.8284D3=4.2426D4=5.6568⇒u0=0/5.6568=0.0000u1=1.4142/5.6568=0.2500u2=2.8284/5.6568=0.5000u3=4.2426/5.6568=0.7500u4=1.0000
步骤 3:计算 B 样条基函数矩阵 AAA
我们将构建矩阵 A∈R5×5A \in \mathbb{R}^{5 \times 5}A∈R5×5,其中每行是对应一个uju_juj,每列是 Ni,k(uj)N_{i,k}(u_j)Ni,k(uj)
使用 Cox-de Boor 公式计算(这里只展示部分关键项)
举例:计算 A2,2=N2,3(u2=0.5)A_{2,2} = N_{2,3}(u_2=0.5)A2,2=N2,3(u2=0.5)
我们知道:
- 节点向量t=[0,0,0,0,0.5,1,1,1,1]t = [0,0,0,0,0.5,1,1,1,1]t=[0,0,0,0,0.5,1,1,1,1]
- 所需节点:t2=0t_2 = 0t2=0, t3=0t_3 = 0t3=0, t4=0.5t_4 = 0.5t4=0.5, t5=1t_5 = 1t5=1, t6=1t_6 = 1t6=1
step1:计算 N2,0(0.5)N_{2,0}(0.5)N2,0(0.5)
t2=0,t3=0⇒t2=t3⇒N2,0(u=0.5)=0t_2 = 0,\quad t_3 = 0 \Rightarrow t_2 = t_3 \Rightarrow N_{2,0}(u=0.5) = 0t2=0,t3=0⇒t2=t3⇒N2,0(u=0.5)=0
因为 t2=t3=0t_2 = t_3 = 0t2=t3=0,没有区间满足 ti≤u<ti+1t_i \le u < t_{i+1}ti≤u<ti+1
step2:计算 N2,1(0.5)N_{2,1}(0.5)N2,1(0.5)
N2,1(u)=u−t2t2+1−t2⋅N2,0(u)+t2+2−ut2+2−t2+1⋅N3,0(u)N_{2,1}(u) = \frac{u - t_2}{t_{2+1} - t_2} \cdot N_{2,0}(u) + \frac{t_{2+2} - u}{t_{2+2} - t_{2+1}} \cdot N_{3,0}(u)N2,1(u)=t2+1−t2u−t2⋅N2,0(u)+t2+2−t2+1t2+2−u⋅N3,0(u)
数值代入:
t2=0,t3=0,t4=0.5t_2 = 0,\quad t_3 = 0,\quad t_4 = 0.5t2=0,t3=0,t4=0.5,
N3,0(0.5)={1,t3≤u<t4⇒0≤0.5<0.50,otherwiseN_{3,0}(0.5) = \begin{cases} 1, & t_3 \le u < t_4 \Rightarrow 0 \le 0.5 < 0.5 \\ 0, & \text{otherwise} \end{cases}N3,0(0.5)={1,0,t3≤u<t4⇒0≤0.5<0.5otherwise
N2,1(0.5)=0+0.5−00.5−0⋅0=0N_{2,1}(0.5) = 0 + \frac{0.5 - 0}{0.5 - 0} \cdot 0 = 0N2,1(0.5)=0+0.5−00.5−0⋅0=0
step3:计算 N2,2(0.5)N_{2,2}(0.5)N2,2(0.5)
递推公式:
N2,2(u)=u−t2t2+2−t2⋅N2,1(u)+t2+3−ut2+3−t2+1⋅N3,1(u)N_{2,2}(u) = \frac{u - t_2}{t_{2+2} - t_2} \cdot N_{2,1}(u) + \frac{t_{2+3} - u}{t_{2+3} - t_{2+1}} \cdot N_{3,1}(u)N2,2(u)=t2+2−t2u−t2⋅N2,1(u)+t2+3−t2+1t2+3−u⋅N3,1(u)
已知 N2,1(u)=0N_{2,1}(u) = 0N2,1(u)=0,计算 N3,1(u)N_{3,1}(u)N3,1(u)
step4:计算 N3,1(0.5)N_{3,1}(0.5)N3,1(0.5)
节点:
- t3=0t_3 = 0t3=0, t4=0.5t_4 = 0.5t4=0.5, t5=1t_5 = 1t5=1
一级基函数:
N3,0(0.5)={1,0≤0.5<0.50,otherwise=0N_{3,0}(0.5) = \begin{cases} 1, & 0 \le 0.5 < 0.5 \\ 0, & \text{otherwise} \end{cases} = 0N3,0(0.5)={1,0,0≤0.5<0.5otherwise=0
N4,0(0.5)={1,0.5≤0.5<11,yes!=1N_{4,0}(0.5) = \begin{cases} 1, & 0.5 \le 0.5 < 1 \\ 1, & \text{yes!} \end{cases} = 1N4,0(0.5)={1,1,0.5≤0.5<1yes!=1
代入:
N3,1(0.5)=0.5−00.5−0⋅0+1−0.51−0.5⋅1=1N_{3,1}(0.5) = \frac{0.5 - 0}{0.5 - 0} \cdot 0 + \frac{1 - 0.5}{1 - 0.5} \cdot 1 = 1N3,1(0.5)=0.5−00.5−0⋅0+1−0.51−0.5⋅1=1
回到 N2,2(0.5)N_{2,2}(0.5)N2,2(0.5)
N2,2(0.5)=0+1−0.51−0⋅1=0.5N_{2,2}(0.5) = 0 + \frac{1 - 0.5}{1 - 0} \cdot 1 = 0.5N2,2(0.5)=0+1−01−0.5⋅1=0.5
step5:计算 N2,3(0.5)N_{2,3}(0.5)N2,3(0.5)
最后一级递推:
N2,3(0.5)=0.5−t2t5−t2⋅N2,2(0.5)+t6−0.5t6−t3⋅N3,2(0.5)N_{2,3}(0.5) = \frac{0.5 - t_2}{t_5 - t_2} \cdot N_{2,2}(0.5) + \frac{t_6 - 0.5}{t_6 - t_3} \cdot N_{3,2}(0.5)N2,3(0.5)=t5−t20.5−t2⋅N2,2(0.5)+t6−t3t6−0.5⋅N3,2(0.5)
上面我们得到 N2,2(0.5)=0.5N_{2,2}(0.5) = 0.5N2,2(0.5)=0.5
再计算 N3,2(0.5)N_{3,2}(0.5)N3,2(0.5)
步骤 4:最小二乘求解控制点 P\mathbf{P}P
step 0 准备数据
控制点个数(未知):5个⇒n+1=55个 ⇒ n+1=55个⇒n+1=5
路径点个数(已知):5个⇒m+1=55个 ⇒ m+1=55个⇒m+1=5
路径点参数 uuu:
u=[0.0,0.25,0.5,0.75,1.0]u = [0.0,\ 0.25,\ 0.5,\ 0.75,\ 1.0]u=[0.0, 0.25, 0.5, 0.75, 1.0]
路径点坐标:
Qx=[01234],Qy=[01012]Q_x = \begin{bmatrix} 0 \\ 1 \\ 2 \\ 3 \\ 4 \end{bmatrix}, \quad Q_y = \begin{bmatrix} 0 \\ 1 \\ 0 \\ 1 \\ 2 \end{bmatrix}Qx=01234,Qy=01012
节点向量(开放均匀):
t=[0,0,0,0,0.5,1,1,1,1]t = [0,\ 0,\ 0,\ 0,\ 0.5,\ 1,\ 1,\ 1,\ 1]t=[0, 0, 0, 0, 0.5, 1, 1, 1, 1]
思路:
我们要解如下线性系统(分别对 x 和 y 方向):
A⋅Px=Qx,A⋅Py=QyA \cdot P_x = Q_x,\quad A \cdot P_y = Q_yA⋅Px=Qx,A⋅Py=Qy
其中:
- AAA 是 5×55 \times 55×5的 B 样条基函数矩阵
- PxP_xPx、PyP_yPy 是控制点的横/纵坐标
step 1:计算基函数矩阵 AAA
我们使用前面画出的基函数图 + 公式,在每个 uju_juj 处计算 5 个基函数:
我们先手动估算这些值(近似值):
uju_juj | N0,3N_{0,3}N0,3 | N1,3N_{1,3}N1,3 | N2,3N_{2,3}N2,3 | N3,3N_{3,3}N3,3 | N4,3N_{4,3}N4,3 |
---|---|---|---|---|---|
0.00 | 1.000 | 0.000 | 0.000 | 0.000 | 0.000 |
0.25 | 0.422 | 0.478 | 0.099 | 0.001 | 0.000 |
0.50 | 0.042 | 0.458 | 0.458 | 0.042 | 0.000 |
0.75 | 0.001 | 0.099 | 0.478 | 0.422 | 0.000 |
1.00 | 0.000 | 0.000 | 0.000 | 0.000 | 1.000 |
所以近似基函数矩阵 A≈A \approxA≈:
A=[1.0000.0000.0000.0000.0000.4220.4780.0990.0010.0000.0420.4580.4580.0420.0000.0010.0990.4780.4220.0000.0000.0000.0000.0001.000]A = \begin{bmatrix} 1.000 & 0.000 & 0.000 & 0.000 & 0.000 \\ 0.422 & 0.478 & 0.099 & 0.001 & 0.000 \\ 0.042 & 0.458 & 0.458 & 0.042 & 0.000 \\ 0.001 & 0.099 & 0.478 & 0.422 & 0.000 \\ 0.000 & 0.000 & 0.000 & 0.000 & 1.000 \\ \end{bmatrix}A=1.0000.4220.0420.0010.0000.0000.4780.4580.0990.0000.0000.0990.4580.4780.0000.0000.0010.0420.4220.0000.0000.0000.0000.0001.000
step 2:解线性系统 A⋅Px=QxA \cdot P_x = Q_xA⋅Px=Qx
写成矩阵方程:
[100000.4220.4780.0990.00100.0420.4580.4580.04200.0010.0990.4780.422000001]⋅[Px0Px1Px2Px3Px4]=[01234]\begin{bmatrix} 1 & 0 & 0 & 0 & 0 \\ 0.422 & 0.478 & 0.099 & 0.001 & 0 \\ 0.042 & 0.458 & 0.458 & 0.042 & 0 \\ 0.001 & 0.099 & 0.478 & 0.422 & 0 \\ 0 & 0 & 0 & 0 & 1 \end{bmatrix} \cdot \begin{bmatrix} P_{x0} \\ P_{x1} \\ P_{x2} \\ P_{x3} \\ P_{x4} \end{bmatrix} = \begin{bmatrix} 0 \\ 1 \\ 2 \\ 3 \\ 4 \end{bmatrix}10.4220.0420.001000.4780.4580.099000.0990.4580.478000.0010.0420.422000001⋅Px0Px1Px2Px3Px4=01234
这个矩阵是严格正定的,我们可以用最小二乘解:
Px=(ATA)−1ATQxP_x = (A^T A)^{-1} A^T Q_xPx=(ATA)−1ATQx
同理得 PyP_yPy
step 3:数值求解(用代码计算)
我将用 NumPy 解出控制点坐标 Px,PyP_x, P_yPx,Py。接下来我来计算。
我们成功求得了 B 样条的控制点坐标(用于拟合路径点 Qj\mathbf{Q}_jQj):
控制点坐标(结果)
横坐标 PxP_xPx:
Px=[−7.2×10−161.58312.41694.00004.0000]P_x = \begin{bmatrix} -7.2 \times 10^{-16} \\ 1.5831 \\ 2.4169 \\ 4.0000 \\ 4.0000 \end{bmatrix}Px=−7.2×10−161.58312.41694.00004.0000
注意:
- 第一个值约为 0(数值误差)
- 最后一个控制点强制等于路径终点 x=4x=4x=4,符合开放均匀节点向量的端点逼近性
纵坐标 PyP_yPy:
Py=[−5.6×10−172.7537−3.24935.40412.0000]P_y = \begin{bmatrix} -5.6 \times 10^{-17} \\ 2.7537 \\ -3.2493 \\ 5.4041 \\ 2.0000 \end{bmatrix}Py=−5.6×10−172.7537−3.24935.40412.0000
可见:
- 控制点纵坐标可能超出了原路径点的范围(例如存在负值和大于2的点)
- 这是最小二乘拟合的结果,为追求整体平滑度而允许个别偏离
步骤 5:构造曲线并采样
在 u∈[0,1]u \in [0,1]u∈[0,1]上均匀采样 100 个点,计算:
C(u)=∑i=0nNi,k(u)⋅Pi\mathbf{C}(u) = \sum_{i=0}^{n} N_{i,k}(u) \cdot \mathbf{P}_iC(u)=∑i=0nNi,k(u)⋅Pi
每一维分别计算:
x(u)=∑i=0nNi,k(u)⋅Px,iy(u)=∑i=0nNi,k(u)⋅Py,ix(u) = \sum_{i=0}^{n} N_{i,k}(u) \cdot P_{x,i} \quad y(u) = \sum_{i=0}^{n} N_{i,k}(u) \cdot P_{y,i}x(u)=∑i=0nNi,k(u)⋅Px,iy(u)=∑i=0nNi,k(u)⋅Py,i
7.为什么 B 样条的“构造 + 采样” = 路径平滑?
B 样条的设计初衷就是“构造一条平滑、连续、可调控的曲线”,它的几大数学特性正好符合路径平滑的核心目标:
路径平滑目标 | B样条如何满足 |
---|---|
避免锯齿(C⁰ 断点) | B 样条曲线是 Ck−1C^{k-1}Ck−1 连续的,三次即 C2C^2C2(有连续曲率) |
避免剧烈弯曲 | 曲线由低阶多项式拼接构成,天然避免“龙格现象” |
局部调整 | 修改一个控制点,只影响其对应的 k+1k+1k+1 个曲线段(局部支持性) |
连续可导与可积性 | 可求导、可计算曲率 κ(u)\kappa(u)κ(u),可作为优化目标 |
数值稳定性 | 每段曲线只依赖局部节点向量和控制点,不易受全局抖动影响 |
7.1 构造形式的回顾
回顾 B 样条曲线的表达式:
C(u)=∑i=0nNi,k(u)⋅Pi\mathbf{C}(u) = \sum_{i=0}^{n} N_{i,k}(u) \cdot \mathbf{P}_iC(u)=∑i=0nNi,k(u)⋅Pi
其中:
- Pi\mathbf{P}_iPi:控制点,决定曲线的形状;
- Ni,k(u)N_{i,k}(u)Ni,k(u):第 iii 个 k阶 B 样条基函数,由节点向量{ti}\{t_i\}{ti} 决定其“形状”和“支持区间”;
- u∈[0,1]u \in [0,1]u∈[0,1]:参数域;
7.2 为什么这样可以实现平滑路径优化?
7.2.1. 基函数提供“天然平滑基底”
B 样条基函数是:
- 分段多项式函数(每段通常为三次);
- 具有连续导数(k次B样条具有 Ck−1C^{k-1}Ck−1连续性);
- 局部支持性:每个 Ni,k(u)N_{i,k}(u)Ni,k(u) 仅在有限区间内非零;
- 满足单位划分:∑iNi,k(u)=1\sum_i N_{i,k}(u) = 1∑iNi,k(u)=1;
这意味着:
- 构造出的曲线天然光滑,不会突变;
- 任一控制点 Pi\mathbf{P}_iPi 只影响有限区间;
- 不易出现高阶多项式拟合常见的振荡问题(如龙格现象);
7.2.2. 控制点的解来自于“最小二乘优化问题”
目标是找到最合适的控制点,使曲线尽可能接近路径点 {Qj}\{Q_j\}{Qj}:
minP0,…,Pn∑j=0m∥C(uj)−Qj∥2\min_{\mathbf{P}_0, \dots, \mathbf{P}_n} \sum_{j=0}^m \| \mathbf{C}(u_j) - \mathbf{Q}_j \|^2minP0,…,Pn∑j=0m∥C(uj)−Qj∥2
这个目标就是:
- 用一个 B 样条构造出的平滑函数,最小化拟合原始离散路径的残差;
- 类似线性回归、傅里叶逼近,只不过这里使用的是 B 样条基底而不是正弦或多项式;
所以,控制点不是随便设的,而是“最优逼近意义下”的最优参数解。
7.2.3. 节点向量控制曲线的局部性和形状
节点向量 {ti}\{t_i\}{ti} 定义了每段曲线的“范围”和“影响因子”,可以设计为:
- 均匀节点:等距平滑;
- 弦长参数化节点:根据原始点的距离差控制;
- 非均匀节点:精细调整,靠近关键拐点提高分辨率;
📌 在优化中,节点向量相当于决定了“每段曲线在哪一段起作用”,而控制点决定了“怎么作用”。
7.3 将这套方法理解为“正则化的逼近问题”
可以从另一个角度理解 B 样条的路径拟合是如下目标:
minC(u)∑j∥C(uj)−Qj∥2+λ∫01∣d2C(u)du2∣2duC(u)\min_{\mathbf{C}(u)} \sum_j \|\mathbf{C}(u_j) - Q_j\|^2 + \lambda \int_0^1 \left|\frac{d^2 \mathbf{C}(u)}{du^2}\right|^2 duC(u)minC(u)∑j∥C(uj)−Qj∥2+λ∫01du2d2C(u)2duC(u)
- 第一项:拟合误差
- 第二项:平滑项,惩罚曲率剧烈变化(即,保持二阶导数小)
- B样条构造出的曲线天然具有小的二阶导数 ⇒ 在无须显式正则化的情况下就具备平滑效果
小结对照
构造元件 | 在路径平滑中起的作用 |
---|---|
控制点Pi\mathbf{P}_iPi | 决定曲线的几何形状,通过最小残差确定 |
节点向量{ti}\{t_i\}{ti} | 决定每个控制点的影响范围和曲线的分段结构 |
基函数 Ni,k(u)N_{i,k}(u)Ni,k(u) | 提供平滑的局部基底,拼接成整体光滑路径 |
8.python实现代码
import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import BSpline# 控制点设置(可以根据需要修改)
control_points = np.array([[0.0, 0.0],[1.0, 2.0],[2.0, -1.0],[3.0, 2.0],[4.0, 0.0],[5.0, 2.0],[6.0, 0.0]
])# 曲线度数 k=3(三次B样条)
k = 3
n = len(control_points) - 1
m = n + k + 1# 构造钳位型均匀节点向量 t(首尾重复k次)
t = np.concatenate((np.zeros(k),np.linspace(0, 1, m - 2 * k + 1),np.ones(k)
))# 提取控制点 x 和 y 坐标
cx = control_points[:, 0]
cy = control_points[:, 1]# 分别构造 x 和 y 的 B 样条函数
spl_x = BSpline(t, cx, k)
spl_y = BSpline(t, cy, k)# 采样参数 u(均匀采样 200 个点)
u = np.linspace(0, 1, 200)# 生成 B 样条曲线坐标
curve = np.vstack((spl_x(u), spl_y(u))).T# 绘图
plt.figure(figsize=(8, 5))
plt.plot(curve[:, 0], curve[:, 1], label='B-spline Curve', lw=2)
plt.plot(control_points[:, 0], control_points[:, 1], 'o--', label='Control Points')
plt.title('B-spline Curve Example')
plt.xlabel('x')
plt.ylabel('y')
plt.axis('equal')
plt.grid(True)
plt.legend()
plt.show()