使用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})")