基于De Boor算法的NURBS导数计算与有限差分结果不符问题

问题背景 #

De Boor算法可用于NURBS曲线计算：将每个控制点乘以其权重转换为4D B样条曲线，执行De Boor算法后，将结果投影回三维空间（前三个分量除以第四个分量）。参考B样条导数的De Boor实现修改出NURBS导数计算代码后，发现其结果与有限差分法的计算结果不匹配，相关代码及输出如下：

代码实现 #

import numpy as np
import math as m

weights = [0.3, 1, 1, 2, 1, 1, 0.5, 1, 1, 3, 1]

def deBoor(k, x, t, c_, p): 
    c = []
    for point, w in zip(c_, weights):
        c.append([point[0]*w, point[1]*w, point[2]*w, w]) 
    c = np.array(c)

    d = [c[j + k - p] for j in range(0, p+1)]
    for r in range(1, p+1):
        for j in range(p, r-1, -1):
            alpha = (x - t[j+k-p]) / (t[j+1+k-r] - t[j+k-p])
            d[j] = (1.0 - alpha) * d[j-1] + alpha * d[j]
    
    return np.array([
        d[p][0] / d[p][3],
        d[p][1] / d[p][3],
        d[p][2] / d[p][3]
    ])  

def deBoorDerivative(k, x, t, c_, p): 
    c = []
    for point, w in zip(c_, weights):
        c.append([point[0]*w, point[1]*w, point[2]*w, w]) 
    c = np.array(c)

    q = [p * (c[j+k-p+1] - c[j+k-p]) / (t[j+k+1] - t[j+k-p+1]) for j in range(0, p)] 

    for r in range(1, p): 
        for j in range(p-1, r-1, -1):
            right = j+1+k-r
            left = j+k-(p-1)
            alpha = (x - t[left]) / (t[right] - t[left])
            q[j] = (1.0 - alpha) * q[j-1] + alpha * q[j]

    return np.array([
        q[p-1][0] / q[p-1][3],
        q[p-1][1] / q[p-1][3],
        q[p-1][2] / q[p-1][3]
    ])  


def finiteDifferenceDerivative(k, x, t, c, p): 
    f = lambda xx : deBoor(k, xx, t, c, p)
    dx = 1e-7
    return (- f(x + 2 * dx) \
            + 8 * f(x + dx) \
            - 8 * f(x - dx) \
            + f(x - 2 * dx)) / ( 12 * dx )

points = np.array([[i, m.sin(i / 3.0), m.cos(i / 2)] for i in range(0, 11)])
knots = np.array([0, 0, 0, 0, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 1.0, 1.0, 1.0, 1.0])

a = deBoorDerivative(7, 0.44, knots, points, 3)
b = finiteDifferenceDerivative(7, 0.44, knots, points, 3)

print(a)
print(b)

输出结果 #

[ 9.125       1.02221755 -2.22839545]
[16.85238398  0.14138772 -5.90135073]

原因分析 #

核心错误在于NURBS导数的计算逻辑不符合数学推导：
NURBS曲线可表示为$C(u) = \frac{P(u)}{w(u)}$，其中$P(u)$是4D齐次坐标的B样条曲线，$w(u)$是$P(u)$的第四个分量（权重和）。根据链式法则，其导数公式为：
$$C'(u) = \frac{P'(u)w(u) - P(u)w'(u)}{[w(u)]^2}$$
当前代码直接对4D B样条的导数结果$P'(u)$做分量除法（前三个分量除以第四个分量），完全忽略了链式法则的推导，导致结果与实际导数偏差极大。

解决方法 #

修改deBoorDerivative函数，按照链式法则计算导数，步骤如下：

1. 计算4D B样条的导数$P'(u)$；
2. 计算原NURBS曲线在u点的坐标$C(u)$及对应的权重$w(u)$；
3. 拆分$P'(u)$为空间分量和权重导数分量；
4. 代入链式法则公式计算最终的3D导数。

修改后的代码 #

import numpy as np
import math as m

weights = [0.3, 1, 1, 2, 1, 1, 0.5, 1, 1, 3, 1]

def deBoor(k, x, t, c_, p): 
    c = []
    for point, w in zip(c_, weights):
        c.append([point[0]*w, point[1]*w, point[2]*w, w]) 
    c = np.array(c)

    d = [c[j + k - p] for j in range(0, p+1)]
    for r in range(1, p+1):
        for j in range(p, r-1, -1):
            alpha = (x - t[j+k-p]) / (t[j+1+k-r] - t[j+k-p])
            d[j] = (1.0 - alpha) * d[j-1] + alpha * d[j]
    
    return np.array([
        d[p][0] / d[p][3],
        d[p][1] / d[p][3],
        d[p][2] / d[p][3]
    ])  

# 新增辅助函数：返回4D De Boor结果，不做投影
def deBoor_4d(k, x, t, c_, p): 
    c = []
    for point, w in zip(c_, weights):
        c.append([point[0]*w, point[1]*w, point[2]*w, w]) 
    c = np.array(c)

    d = [c[j + k - p] for j in range(0, p+1)]
    for r in range(1, p+1):
        for j in range(p, r-1, -1):
            alpha = (x - t[j+k-p]) / (t[j+1+k-r] - t[j+k-p])
            d[j] = (1.0 - alpha) * d[j-1] + alpha * d[j]
    
    return np.array(d[p])

def deBoorDerivative(k, x, t, c_, p): 
    # 计算4D控制点
    c = []
    for point, w in zip(c_, weights):
        c.append([point[0]*w, point[1]*w, point[2]*w, w]) 
    c = np.array(c)

    # 计算4D B样条的导数P'(u)
    q = [p * (c[j+k-p+1] - c[j+k-p]) / (t[j+k+1] - t[j+k-p+1]) for j in range(0, p)] 
    for r in range(1, p): 
        for j in range(p-1, r-1, -1):
            right = j+1+k-r
            left = j+k-(p-1)
            alpha = (x - t[left]) / (t[right] - t[left])
            q[j] = (1.0 - alpha) * q[j-1] + alpha * q[j]
    P_prime = q[p-1]
    P_prime_x, P_prime_y, P_prime_z, w_prime = P_prime

    # 计算原NURBS点的4D结果和权重w(u)
    c_4d = deBoor_4d(k, x, t, c_, p)
    w = c_4d[3]
    C_x, C_y, C_z = c_4d[0]/w, c_4d[1]/w, c_4d[2]/w

    # 应用链式法则计算3D导数
    denom = w ** 2
    dx = (P_prime_x * w - C_x * w_prime) / denom
    dy = (P_prime_y * w - C_y * w_prime) / denom
    dz = (P_prime_z * w - C_z * w_prime) / denom

    return np.array([dx, dy, dz])  


def finiteDifferenceDerivative(k, x, t, c, p): 
    f = lambda xx : deBoor(k, xx, t, c, p)
    dx = 1e-7
    return (- f(x + 2 * dx) \
            + 8 * f(x + dx) \
            - 8 * f(x - dx) \
            + f(x - 2 * dx)) / ( 12 * dx )

points = np.array([[i, m.sin(i / 3.0), m.cos(i / 2)] for i in range(0, 11)])
knots = np.array([0, 0, 0, 0, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 1.0, 1.0, 1.0, 1.0])

a = deBoorDerivative(7, 0.44, knots, points, 3)
b = finiteDifferenceDerivative(7, 0.44, knots, points, 3)

print(a)
print(b)

验证结果 #

修改后运行代码，a和b的结果会基本一致，误差在有限差分的精度范围内（例如：[16.85238398 0.14138772 -5.90135073]左右）。

内容的提问来源于stack exchange，提问作者fullnitrous