Search code examples
pythonnumpyspline

Singular matrix during B Spline interpolation

According to the literature about B Splines, including Wolfram Mathworld, the condition for Cox de Boor's recursive function states that:

enter image description here

In python, this would translate to:

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

where:

  • d_: degree of the curve
  • knots_: knot vector
  • k_: index of the knot
  • t_: parameter value {0.0,...,1.0} (reparametrized)

However, this seems to generate a Singular matrix, when creating the linear system intended for interpolation, not approximation. For example, with 4 points:

  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 I change the second part of the condition to:

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

I get the correct matrix and the correct spline.

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

Why does the code need to deviate from the mathematical condition in order to get the correct results and the iterator reaching the last element?

See below the complete example code:

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 )

Result:

enter image description here

Solution

  • The mathematical condition is the correct one, B-Spline basis functions are defined on a half-open interval.

    However, it presents a problem when the knot vector is clamped, as you show in your example. The problem is that in this case the mathematical function isn't defined at t=1, and the result in the code is that it evaluates to 0 as you showed.

    What you would like to have in this case, is for the function to be evaluated to 1 (which is the limit at t=1). One way to achieve this is to apply your fix of using <= instead of <.

    The figures below show a plot of the resulting basis function for n=5 (for which your default knot vector is [0,0,0,0,0.5,1,1,1,1]) using your fix and using the original mathematical definition. The functions were sampled at t = linspace(0,1,101).

    enter image description here enter image description here

    As can be seen, both implementations result in similar basis functions in "almost all" values.

    However, while the <= fix sets the value at t=1 to 1 as we wanted, it clearly ruins the functions at the inner knot value t=0.5. The original implementation, on the other hand, ruins the function at the end knot value t=1 as we already know, but only there and only for the last function.

    Basically, one can say that the basis functions implemented with the <=-fix are "almost-correct" except for inner-knot values. Here is another example figure to demonstrate the phenomeneon for n=7.

    enter image description here

    So, now that we know what the problem is and what we want the result to be, one can think of different ways to fix it. The following code is an example of such a fix, it wraps the recursive definition with a function that checks for (d+1)-multiplicity knots and sets it to 1 (a clamped knot vector is the most common use-case of this).

    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

    And here is the result we get when using this function on the example from above.

    enter image description here