奇数点插值样条异常

Question

我实现了三次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个点（控制顶点）时结果OK：

插值 5 个点（控制顶点）时结果不正常：

插值6个点（控制顶点）时结果OK：

等等...

我注意到两件事：

1. 当控制顶点数为奇数时，样条线无法正确插值。
2. 当度数变为二次方时，样条曲线可以对任意数量的顶点进行正确插值。因此，如果您更改

   d = 2

   ，在

   B_spline

   函数中，曲线将为奇数和偶数控制顶点正确插值。

cox de boor函数是正确的，并且根据数学表达式，但对第二个条件表达式

t[i] <= t **<=** t[i+1]

进行了小改动（请参阅我之前的SO问题here了解更多详细信息）。另外，我使用 numpy 来求解线性系统，它也按预期工作。除了

np.linalg.solve

之外，我还尝试过

np.linalg.lstsq

，但它返回相同的结果。

老实说，我不知道该将这种异常行为归咎于何处。什么可能导致此问题？

Answer 1

所描述的异常行为非常有趣，其原因也很微妙。 基本上，此行为的根本原因是 Cox-DeBoor 函数实现和

<=

修复。

在我对OP之前的SO问题的回答中，我详细解释了为什么这个修复是错误的。 简而言之，该实现构造了除内结值之外“几乎正确”的基函数。 在内结之外的间隔处，基函数给出正确的值 - 包括 t=1

 值，
这就是首先进行修复的原因（以防止插值矩阵中出现零行）。

有了这些知识，我们就可以解释这种神秘的异常行为是如何发生的。

需要注意的主要事情是数组

ts

 和 

knots

，它们作为参数传递给 

interpolate()

 函数
在函数 

B_spline()

 内。

interpolate()

 函数评估所有 

ts

 的基函数。
如果 

ts

 的任何内部值都不等于内部结，那么所有这些值都将是正确的，结果也将是正确的。
然而，如果存在任何

ts[i]==inner_knot

，那么那里的值就会错误并会破坏结果。

现在，剩下的就是尝试解释问题中指出的两件事。

首先，请注意结已夹紧且间距相等 - 这就是

knotvector.generate( d, n )

 的作用。
因此，对于 

d=3

，结是：

[0,0,0,0, 1/(n-3), 2/(n-3),..., (n-4)/(n-3), 1,1,1,1]

。
例如，对于 

n=5

，结向量为 

[0,0,0,0,0.5,1,1,1,1]

。
更一般地，可以使用命令 

knots[d:-d] = np.linspace(0, 1, n-d+1)

 计算内结（包括 0 和 1）。

ts

 使用命令 

np.linspace( 0.0, 1.0, n )

（即 

[0, 1/(n-1), 2/(n-1),..., (n-2)/(n-1), 1]

）计算。

“当控制顶点数为奇数时，样条曲线无法正确插值”：

这是因为当

n

 是奇数并且 

d=3

 时，

n-d+1

 也是奇数，那么 0.5 一定是结之一。
按照同样的道理，0.5 也一定是 

ts

 之一。
因此，我们得到了错误的答案！

有趣的是，如果

n

是

偶数和d=3

，就可以证明

n-3

和

n-1

之间不存在公因数
（

两个连续的奇数没有公因数） 因此 ts

 的内部值不等于内部结，结果将是正确的！
例如，对于 

n=6

、

ts = [0, 1/5, 2/5, 3/5, 4/5, 1]

 和 

knots = [..,0, 1/3, 2/3, 1,..]

，所以没有
骨折 

i/3

 等于骨折 

j/5

。

“当度数变为二次时，样条曲线可以对任意数量的顶点进行正确插值”：

这是因为当

d=2

时，内结将是

np.linspace(0, 1, n-1)

（即

[0,0,0, 1/(n-2), 2/(n-2),..., (n-3)/(n-2), 1,1,1]

）并且

ts

将是

np.linspace(0, 1, n)

。
由于

n-2

和

n-1

之间没有公因数，所以我们不会有

i/(n-2) == j/(n-1)

，所以总是

ts != knots

，结果将是正确的！

我相信这解释了这种有趣的异常插值行为。