B样条曲线在给定区间内曲率极值的计算方法

1. B样条曲线定义

设 p 次 B 样条曲线为：

C(u) = Σᵢ₌₀ⁿ Pᵢ Nᵢ,ₚ(u), u ∈ [uₐ, u_b] ⊆ [uₚ, uₘ₋ₚ]

其中：

• Pᵢ ∈ ℝᵈ（d = 2 或 3）为控制点；
• Nᵢ,ₚ(u) 为由节点向量 U = {u₀, u₁, …, uₘ} 定义的 p 次 B 样条基函数。

2. 曲率表达式

（1）平面曲线（d = 2）

曲率 κ(u) 为：

κ(u) = |C′(u) × C″(u)| / ‖C′(u)‖³

其中：

• C′(u) 和 C″(u) 分别为一阶和二阶导数；
• “×” 表示二维向量的伪叉积：若 C′ = (x′, y′)，C″ = (x″, y″)，则
  C′ × C″ = x′ y″ − y′ x″（标量）；
• ‖·‖ 为欧几里得范数。

（2）空间曲线（d = 3）

曲率 κ(u) 为：

κ(u) = ‖C′(u) × C″(u)‖ / ‖C′(u)‖³

其中 “×” 为三维向量叉积，“‖·‖” 为向量模长。

注：本文以平面情形为主，空间情形类似，仅将叉积结果取模。

3. 导数计算（B样条导数性质）

B 样条曲线的 k 阶导数仍为 B 样条形式。一阶和二阶导数可表示为：

C′(u) = p Σᵢ₌₀ⁿ⁻¹ (Pᵢ₊₁ − Pᵢ) / (uᵢ₊ₚ₊₁ − uᵢ₊₁) · Nᵢ,ₚ₋₁(u)

C″(u) = p(p−1) Σᵢ₌₀ⁿ⁻² Qᵢ · Nᵢ,ₚ₋₂(u)

其中 Qᵢ 为二阶差商控制点（可通过递推或自动微分获得）。
实际计算中，可使用 de Boor 算法或其导数扩展高效求 C′(u) 和 C″(u)。

4. 曲率极值条件

曲率极值出现在以下两类点中：

1. 驻点：满足 dκ/du = 0 的内点 u ∈ (uₐ, u_b)；
2. 端点：u = uₐ 或 u = u_b。

因此，需解方程：

dκ(u)/du = 0

由于 κ(u) 是有理函数（分子为多项式，分母为范数的三次方），其导数为高阶有理函数，一般无解析解，需数值求解。

5. 数值求解步骤

步骤 1：参数区间离散化
在 [uₐ, u_b] 内生成初始采样点，如均匀网格或基于曲率变化自适应采样。

步骤 2：计算曲率及其导数（或差分）
对每个 u，用 B 样条求值算法计算：

• C′(u), C″(u)
• κ(u) = |x′ y″ − y′ x″| / (x′² + y′²)^{3/2}

步骤 3：定位候选极值点

• 检查 κ(u) 在采样点中的局部极大/极小；
• 对每个疑似极值区间 [uᵢ, uᵢ₊₁]，使用数值优化方法（如 Brent 法、牛顿法）求解 dκ/du = 0。
• 实用技巧：可最小化/最大化目标函数 f(u) = κ(u)，使用无导数优化器（如 golden-section search）避免显式求导。

步骤 4：比较端点与内部极值
最终曲率极值为：

• κ_max = max{ κ(uₐ), κ(u_b), κ(u₁*), κ(u₂*), … }
• κ_min = min{ κ(uₐ), κ(u_b), κ(u₁*), κ(u₂*), … }

其中 uⱼ* 为内部驻点。

6. 注意事项

• 当 ‖C′(u)‖ ≈ 0 时（曲线速度接近零），曲率可能发散，需排除奇异点；
• 若 B 样条为非均匀有理 B 样条（NURBS），需先转换为有理形式或使用齐次坐标求导；
• 对于高精度需求，建议使用区间分析或贝塞尔 clipping 方法保证极值不被遗漏。