使用deboor法计算三次B样条曲线在参数为u处的位置的方法介绍

核心概念

1. 控制点 (Control Points): P[0], P[1], ..., P[n]
2. 节点向量 (Knot Vector): U = [u_0, u_1, ..., u_{m}]。对于定义在区间 [a, b] 上的 p 次B样条（这里是3次），通常需要 m+1 = n + p + 2 个节点。
   • 常见的均匀周期性节点向量（假设区间 [0, 1]）：U = [0, 0, 0, 0, 1/(n-2), 2/(n-2), ..., 1, 1, 1, 1]。注意首尾各有 p+1=4 个重复节点。
   • 或者更通用的非周期性均匀节点向量：U = [0, 0, 0, 0, 1, 2, ..., n-2, n-2, n-2, n-2] / (n-2) (如果 n > 2)。
3. 次数 (Degree): p = 3 (三次)。
4. De Boor 算法: 通过递归地线性插值控制点来计算 C(u)。

De Boor 算法步骤

目标是计算 C(u)，即曲线上参数为 u 的点。

1. 找到正确的区间 (Find Span): 找到节点向量中满足 u_k <= u < u_{k+1} 的索引 k。这个 k 就是“节点跨度”(knot span)。由于首尾有重复节点，实际有效的 k 范围是 [p, m-p-1] 或 [3, len(U)-4]。

   • 特殊情况: 如果 u 等于最后一个节点 u_m，则令 k = m - p - 1。

2. 初始化控制点 (Initialize Control Points): 取与节点跨度 k 相关的 p+1=4 个控制点作为初始的 De Boor 点： d_i^(0) = P_{k-p+i}，其中 i = 0, 1, 2, 3。 即：d0 = P[k-3], d1 = P[k-2], d2 = P[k-1], d3 = P[k]。

3. 递归计算 (Recursive Calculation): 对于 r = 1 到 p (即 r=1, 2, 3)，以及对于 i = p-r 到 p (随着 r 增加，i 的范围缩小)，计算： d_i^(r) = alpha_i^(r) * d_{i-1}^(r-1) + (1 - alpha_i^(r)) * d_i^(r-1) 其中插值系数 alpha_i^(r) 为： alpha_i^(r) = (u - u_{k-p+i}) / (u_{k+i+1-r} - u_{k-p+i})

   这个过程会将4个点逐步缩减为1个点。

4. 结果: 经过 p=3 轮迭代后，剩下的唯一一个点 d_p^(p) = d_3^(3) 就是 C(u)。

Python代码实现：

import numpy as np

def find_span(n, p, u, U):
"""
找到参数 u 在节点向量 U 中的跨度 k。

参数:
n: 控制点数量减一 (即控制点索引最大值)
p: B样条次数
u: 参数值
U: 节点向量 (长度 m+1 = n+p+2)

返回:
k: 节点跨度索引
"""
# 使用二分查找确保效率
if u >= U[n+1]:
return n

low = p
high = n + 1
mid = (low + high) // 2

while u < U[mid] or u >= U[mid + 1]:
if u < U[mid]:
high = mid
else:
low = mid
mid = (low + high) // 2
return mid

def de_boor(k, x, d, i, p, U):
"""
递归计算De Boor点。这里用循环实现更高效。
"""
# 不采用递归，而是用嵌套循环
pass # 此函数不直接使用，逻辑在主函数中

def evaluate_bspline_point(P, U, u, p=3):
"""
使用De Boor算法计算三次B样条曲线在参数u处的点。

参数:
P: 控制点列表，每个点可以是元组或列表，如 [(x0,y0), (x1,y1), ...]
U: 节点向量 (list or numpy array)
u: 参数值
p: B样条次数 (默认为3)

返回:
C_u: 曲线上对应u的点 (numpy array)
"""
n = len(P) - 1  # 控制点索引最大值
m = len(U) - 1  # 节点向量索引最大值

# 检查次数和节点向量长度
assert m == n + p + 1, f"节点向量长度错误。期望 m={n+p+1}, 实际 m={m}"

# 1. 找到跨度 k
k = find_span(n, p, u, U)

# 2. 初始化De Boor点 (只取相关的 p+1 个点)
# d 是一个 (p+1) x dim 的数组，用于存储每轮迭代的结果
dim = len(P[0])  # 点的维度 (2D, 3D等)
d = np.zeros((p+1, dim))

for j in range(p+1):
d[j] = np.array(P[k - p + j])

# 3. 递归计算 (使用循环)
for r in range(1, p+1):  # r = 1, 2, 3
# 当前轮次需要计算的点的索引范围: 从 (p-r) 到 p
for i in range(p, r-1, -1): # i 从 p 递减到 r
# 计算 alpha
left_idx = k - p + i
right_idx = k + i + 1 - r
denominator = U[right_idx] - U[left_idx]

if abs(denominator) < 1e-12:  # 避免除以零
alpha = 0.0
else:
alpha = (u - U[left_idx]) / denominator

# 线性插值
d[i] = alpha * d[i-1] + (1.0 - alpha) * d[i]

# 4. 结果在 d[p] 中
return d[p]

# --- 示例 ---
if __name__ == "__main__":
# 定义控制点 (5个点，n=4)
control_points = [
(0, 0),
(1, 2),
(3, 1),
(4, 4),
(6, 0)
]

# 构造均匀非周期性节点向量 (p=3, n=4 -> m = 4+3+1=8, U长度9)
# U = [0,0,0,0, 1, 2, 3, 3,3,3] 归一化到 [0,1] 区间
# 这里为了简单，直接构造
# 有效内部节点数: n - p = 1, 所以只有一个内部节点
# 常见做法: U = [0,0,0,0, 0.5, 1,1,1,1] 或 [0,0,0,0,1,1,1,1] 但后者退化
# 我们用: U = [0,0,0,0, 1, 2, 3, 3,3,3] / 3 -> [0,0,0,0, 1/3, 2/3, 1,1,1,1]
# 但 n=4, p=3, m=8, 所以需要9个节点
# 更标准: U = [0,0,0,0, 1, 2, 2,2,2] 如果区间是[0,2]，但通常归一化
# 简单起见，使用对称的: 
# U = [0,0,0,0, 1, 1,1,1] 但这是退化的 (所有内部节点相同)
# 正确构造: 对于 n=4, p=3, 均匀非周期: U = [0,0,0,0, 1, 2, 2,2,2] (如果区间[0,2])
# 或者归一化: U = [0,0,0,0, 0.5, 1,1,1,1]

# 使用: U = [0,0,0,0, 0.5, 1,1,1,1]
knot_vector = [0, 0, 0, 0, 0.5, 1, 1, 1, 1]

# 计算几个点
u_values = [0.0, 0.25, 0.5, 0.75, 1.0]

print("控制点:", control_points)
print("节点向量:", knot_vector)
print("\n参数 u 和对应的曲线点:")

for u in u_values:
point = evaluate_bspline_point(control_points, knot_vector, u, p=3)
print(f"u={u:.2f} -> C(u)=({point[0]:.3f}, {point[1]:.3f})")