奇数控制点的三次B样条插值异常问题排查

问题代码 #

import numpy as np
import math
from geomdl import knotvector

def cox_de_boor( d_, t_, k_, knots_):
    if (d_ == 0):
        if ( knots_[k_] <= t_ <= knots_[k_+1]):
            return 1.0
        return 0.0
    
    denom_l = (knots_[k_+d_] - knots_[k_])
    left = 0.0    
    if (denom_l != 0.0):
        left = ((t_ - knots_[k_]) / denom_l) * cox_de_boor(d_-1, t_, k_, knots_)
    
    denom_r = (knots_[k_+d_+1] - knots_[k_+1])
    right = 0.0
    if (denom_r != 0.0):
        right = ((knots_[k_+d_+1] - t_) / denom_r) * cox_de_boor(d_-1, t_, k_+1, knots_)

    return left + right

def interpolate( d_, P_, n_, ts_, knots_ ):
    A = np.zeros((n_, n_))
    
    for i in range(n_):
        for j in range(n_):
            A[i, j] = cox_de_boor(d_, ts_[i], j, knots_)
    
    control_points = np.linalg.solve(A, P_)
    return control_points

def create_B_spline( d_, P_, t_, knots_):
    sum = Vector()   # just a vector class.
    for i in range( len(P_) ):
        sum += P_[i] * cox_de_boor(d_, t_, i, knots_)
    return sum
    
def B_spline( points_ ):
    d     = 3   # change to 2 for quadratic.
    P     = np.array( points_ )
    n     = len( P )
    ts    = np.linspace( 0.0, 1.0, n )
    knots = knotvector.generate( d, n ) # len = n + d + 1
    
    control_points = interpolate( d, P, n, ts, knots)
    
    crv_pnts = []
    for i in range(10):
        t = float(i) / 9
        crv_pnts.append( create_B_spline(d, control_points, t, knots) )
    return crv_pnts


control_points = [ [float(i), math.sin(i), 0.0] for i in range(4) ]
cps = B_spline( control_points )

测试现象 #

• 插值4、6个控制点（插值点）时结果正常
• 插值5个等奇数个控制点时结果异常
• 将样条阶数改为二次（d=2）后，无论控制点数量奇偶，插值结果均正常

已排查内容 #

• Cox-de-Boor函数符合数学定义（仅调整了条件表达式写法）
• 尝试np.linalg.solve和np.linalg.lstsq求解线性系统，结果一致

1. 均匀参数与节点向量的冲突 #

你使用geomdl.knotvector.generate(d, n)生成的是标准均匀节点向量：对于三次B样条（d=3），节点向量两端重复d+1次（即4个0和4个1），中间内部节点均匀分布。当插值点数量n为奇数时，np.linspace(0.0, 1.0, n)生成的均匀插值参数中，中间点恰好落在内部节点的位置。

以n=5为例：

• 节点向量为[0,0,0,0, 0.5, 1,1,1,1]
• 插值参数ts为[0, 0.25, 0.5, 0.75, 1]，其中0.5正好是内部节点

此时在t=0.5处计算三次B样条基函数，会导致插值矩阵A的对应行出现线性相关，矩阵变得奇异（或接近奇异），求解线性方程组时会得到异常的控制点。而偶数个插值点时，ts的取值刚好避开内部节点，矩阵A非奇异，求解正常。

2. 二次B样条无异常的原因 #

当阶数改为二次（d=2）时，即使插值参数点落在内部节点上，二次B样条的基函数在节点处保持C¹连续，不会触发矩阵奇异的情况。同时二次样条的节点向量内部节点间隔与插值参数的匹配关系更稳定，因此无论插值点数量奇偶，求解都能正常进行。

解决方法 #

• 调整插值参数：避免使用严格均匀的linspace，可将中间参数点轻微偏移，或改用弦长参数化（根据插值点间的距离分配参数值），让参数点避开内部节点。
• 使用非均匀节点向量：手动生成节点向量，使内部节点位置与插值参数点错开。

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