B 样条插值期间的奇异矩阵

Question

根据有关 B 样条曲线的文献，包括 Wolfram Mathworld，Cox de Boor 递归函数的条件指出：

在Python中，这将翻译为：

if (d_ == 0):
   if ( knots_[k_] <= t_ < knots_[k_+1]):
      return 1.0
   return 0.0

地点：

• d_

  ：曲线的度数

• knots_

  ：结向量

• k_

  ：结的索引

• t_

  ：参数值{0.0,...,1.0}（重新参数化）

但是，在创建用于“插值”而不是“近似”的线性系统时，这似乎会生成一个奇异矩阵。例如，有 4 点：

  A = [[1.         0.         0.         0.        ]
       [0.2962963  0.44444444 0.22222222 0.03703704]
       [0.03703704 0.22222222 0.44444444 0.2962963 ]
       [0.         0.         0.         **0.**    ]] //The last element (bottom-right) should have been 1.0
  # Error: LinAlgError: file C:\Users\comevo\AppData\Roaming\Python\Python310\site-packages\numpy\linalg\_linalg.py line 104: Singular matrix

如果我将条件的第二部分更改为：

if (d_ == 0):
   if ( knots_[k_] <= t_ <= knots_[k_+1]): // using <= instead of <
      return 1.0
   return 0.0

我得到了正确的矩阵和正确的样条线。

A = [[1.         0.         0.         0.        ]
     [0.2962963  0.44444444 0.22222222 0.03703704]
     [0.03703704 0.22222222 0.44444444 0.2962963 ]
     [0.         0.         0.         1.        ]] // OK

为什么代码需要偏离数学条件才能得到正确的结果并且迭代器到达最后一个元素？

请参阅下面完整的示例代码：

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 = MVector()
    for i in range( len(P_) ):
        sum += P_[i] * cox_de_boor(d_, t_, i, knots_)
    return sum
    
def B_spline( points_ ):
    d     = 3
    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(8) ]
cps = B_spline( control_points )

结果：

数学条件是正确的，B 样条基函数是在

半开

Answer 1

但是，当结向量被夹紧时，就会出现问题，如示例所示。 问题是在这种情况下数学函数没有定义在 t=1， 代码中的结果是它的计算结果为 0，如您所示。

在这种情况下，您希望函数的计算结果为 1（这是

t=1

处的极限）。 实现此目的的一种方法是应用使用

<=

而不是

<

的修复。下图显示了

n=5

的结果基函数图（默认结向量是

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

）使用您的修复并使用原始的数学定义。 函数在

t = linspace(0,1,101)

处采样。

可以看出，两种实现都会在“几乎所有”值中产生类似的基函数。 但是，虽然

<=

修复将

t=1

的值设置为我们想要的 1，但它显然 破坏了内结值

t=0.5

处的功能。 另一方面，正如我们所知，最初的实现破坏了末尾结值

t=1

处的函数， 但仅限于最后一个功能。基本上，可以说使用

<=

修复实现的基函数“几乎正确”，除了内结值之外。 这是另一个示例图来演示

n=7

的现象。

所以，既然我们知道问题是什么以及我们想要的结果是什么，我们就可以想出不同的方法来解决它。 以下代码是此类修复的示例，它使用检查函数来包装递归定义 (d+1)-多重结并将其设置为 1（夹紧结向量是最常见的用例）。

def cox_de_boor(d_, t_, k_, knots_):
    # Handling end-parameter of clamped knot vector (and also the not-so-useful case of (d+1)-multiplicity inner knots)
    if t_ == knots_[k_ + 1] and t_ == knots_[k_ + d_ + 1]:
        return 1.0

    return cox_de_boor_recursive(d_, t_, k_, knots_)


def cox_de_boor_recursive(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_recursive(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_recursive(d_ - 1, t_, k_ + 1, knots_)

    return left + right

这是我们在上面的示例中使用此函数时得到的结果。